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warhead  fragments  from  missiles  killed  at  very  short  range  to  preclude  significant  damage  to  the 
defending  ship.  Furthermore,  the  barrier  would  defeat  the  fuzing  and  structure  of  ASMs  that  have 
penetrated  the  inner  self-defense  layer. 

The  main  thrust  of  this  report  is  to  provide  validations  of  a  hydrodynamics  computer  code 
for  predicting  explosion  plume  behavior.  The  results  of  experiments  conducted  at  NSWC 
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ABSTRACT 


A  mathematical  model  and  computational  code  are  validated  for  predicting  shallow  depth 
explosion  plume  behavior.  The  model  is  based  on  a  generalized  formulation  of  hydrodynamics 
and  uses  an  incompressible  liquid  assumption.  This  formulation  is  well  suited  for  predicting  long¬ 
time  bubble  and  plume  dynamics.  Initial  conditions  for  the  model  are  derived  from  spherically 
symmetric  bubble  theory,  combined  with  empirical  measurements.  The  effects  of  the  spray  dome 
caused  by  the  reflection  of  the  initial  shock  wave  off  the  free  surface  is  modeled  empirically  as  a 
recess  in  the  surface  above  the  charge.  Comparisons  to  photographs  of  experiments  provide 
qualitative  agreement.  Quantitative  measurements  of  plume  heights  and  plume  densities  using 
conductivity  probes  and  microwave  absorption  are  also  compared  to  the  computational  data. 
Results  for  both  single-  and  multiple-charge  shots  are  included. 
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CHAPTER  1 
INTRODUCTION 

During  the  past  50  years  a  large  amount  of  research  has  been  focused  on  the  phenomenon  of 
plumes  generated  by  underwater  explosions.  Reviews  of  much  of  this  work  can  be  found  in  the 
book  by  Cole,1  and  the  reports  of  Young.2,3  These  references  provide  observations  and  empirical 
formulas  for  large  charges,  which  describe  the  general  shape  of  the  resulting  plume,  in  particular,  its 
height  and  width.  Photographs  of  shallow  depth  spark  induced  bubbles  presented  by  Blake  and 
Gibson,4  clearly  show  a  water  jet  rising  above  the  surface  during  the  first  collapse  of  the  bubble. 
Some  interesting  experiments  were  performed  by  Kedrinskii 5  who  suspended  an  inverted  flask  with 
a  “trap”  above  small  charges  and  measured  the  total  amount  of  water  in  the  jets  produced  as  a 
function  of  the  charge  depth.  Besides  this,  very  little  is  known  quantitatively  about  the  total  amount 
water  in  the  plume  as  well  as  the  structure  of  the  plume.  The  lack  of  detailed  plume  knowledge  is 
due  primarily  to  the  opaqueness  of  the  spray  caused  by  the  initial  “spalling”  of  the  surface  from 
the  shock  and  the  breakup  of  the  ejected  plumes.  It  is  also  due  in  part  to  a  lack  of  efficient  and 
robust  computational  procedures  capable  of  predicting  this  complex  phenomenon.  The  purpose  of 
this  report  is  to  begin  to  provide  details  of  the  plume  structure  using  both  experimental  measure¬ 
ments  as  well  as  computational  predictions. 

As  described  by  Colee,1  and  updated  with  a  more  complete  understanding  by  Snay6  and 
Kedrinskii,5  there  are  several  distinct,  but  interrelated  phenomena  typical  to  shallow  depth  explo¬ 
sions.  Upon  detonation  of  the  explosive,  a  spherical  shock  wave  moving  away  from  the  charge  is 
emitted.  This  wave  reflects  off  the  surface  as  a  rarefaction  wave  which  travels  back  down  through 
the  gas  globe  of  detonation  products.  Due  to  the  tension  created  behind  the  rarefaction  wave,  cavita¬ 
tion  can  occur  and  some  of  the  water  can  be  spalled  upwards,  creating  a  “spray  dome.”  It  is  the 
height  and  width  of  this  spray  dome  or  “cupola”  that  is  empirically  fit  to  data  by  Kelsky  et  al.7  and 
Young,2  at  least  at  the  early  times. 

The  greatest  interest  of  this  report  is  the  eruption  of  vertical  and  radial  plumes  caused  by  the 
dynamics  of  the  pulsating  bubble.  This  phase  of  the  problem  has  also  been  controversial  in  conjec¬ 
tures  of  its  cause,  structure,  and  density.  In  particular,  consider  the  case  when  the  scaled  depth 

C  =  — - ,  where  d  is  the  depth  and  Amax  is  the  maximum  bubble  radius,  is  less  than  one. 

^max 

Through  the  use  of  potential  flow  theory,  as  well  as  photographic  evidence  of  spark  bubbles  (in 
which  there  are  little  or  no  shock  effects),  it  has  been  clearly  demonstrated  that  the  initial  or  pri¬ 
mary  plume  consists  of  a  central  column  made  up  almost  entirely  of  water  rising  above  the 
surface.4,5, 6,8  For  large  explosives,  the  velocity  of  the  liquid  jet  rising  above  the  surface  has  been 
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underpredicted  by  potential  flow  procedures.  This  discrepancy  is  very  likely  due  to  the  initial  spalla¬ 
tion  of  the  water  which  leaves  an  indentation  in  the  liquid  surface  above  the  bubble.  Potential  flow 
theory  has  also  shown  that  when  an  indented  surface  is  assumed  for  the  initial  condition  then  the  jet 
is  intensified  in  its  velocity.5 

For  values  of  C  >  0.5,  it  has  been  observed  that  in  the  case  of  spark  generated  bubbles,  the 
free  surface  is  not  broached  by  the  bubble  during  its  first  complete  pulsation.4  Secondary  plumes 
ejected  radially  above  the  surface  after  the  bubble  collapse  have  been  also  been  observed  and  are 
due  to  the  second  expansion  of  the  bubble,  which  at  this  time  forms  an  annular  region.3  It  will  be 
demonstrated  in  this  report  that  these  secondary  plumes  can  be  reliably  simulated  by  our  computa¬ 
tional  method. 

This  report  contains  data  from  small-scale  explosion  tests  performed  at  the  Naval  Surface  War¬ 
fare  Center  -  Carderock  Division  (NSWCCD),  Bethesda,  Maryland,  Test  Pond  in  May  1993,  and  at 
larger  scales  at  Hi-Test  Laboratories  in  a  quarry  in  Arvonia,  Virginia  in  August  1993,  and  at  the 
Briar  Point  facility  of  the  Abdereen  Proving  Grounds,  Aberdeen,  Maryland,  in  June  1994.  The  Car¬ 
derock  tests  were  photographed  with  cameras  both  above  and  below  the  surface  of  the  water,  so  that 
the  motion  of  the  bubble  could  be  observed  in  conjunction  with  the  plume  formation.  Data  from  the 
Arvonia  and  Aberdeen  tests  came  not  only  from  above  surface  photographs,  but  also  from  conduc¬ 
tivity  probes  (held  up  by  cables  above  the  charges)  and  microwave  absorption  measurements 
through  the  plumes.  In  addition  to  using  these  measurements  to  demonstrate  the  feasibility  of  using 
the  plume  as  an  effective  barrier,  they  are  used  to  provide  benchmark  comparisons  for  validating 
numerical  simulations. 

A  description  of  the  numerical  code  used  for  the  plume  predictions  is  given  in  Chapter  2.  The 
validity  of  the  code  on  various  benchmark  problems  is  presented  in  Chapter  3  including  a  com¬ 
parison  of  the  computations  to  the  photographs  of  spark  bubble  experiments  from  Reference  4. 
Some  theory  of  underwater  bubbles  and  the  effects  of  the  free  surface  are  also  given.  This  theory  is 
used  together  with  shallow-depth  period  and  bubble-size  measurements  to  derive  empirical  parame¬ 
ters  for  the  charge  types  under  consideration.  These  new  parameters  are  more  consistent  with  the 
observed  data  and  can  be  expected  to  provide  for  more  accurate  initial  conditions.  Comparisons  to 
the  measured  data  together  with  the  results  of  the  numerical  simulations  are  presented  in  Chapter  4. 
A  brief  summary  of  the  results  and  conclusions  is  contained  in  Chapter  5. 
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CHAPTER  2 

THE  NUMERICAL  METHOD 

The  numerical  simulations  used  in  this  study  were  produced  using  the  code  BUB2D  for  both 
axially  symmetric  and  two-dimensional  problems.  Both  codes  use  an  incompressible  assumption  for 
the  regions  containing  water.  One  advantage  of  such  an  approach  is  that  the  time-step  limitation  is 
based  only  on  the  fluid  velocity  as  compared  to  compressible  formulations  that  have  time-step  limi¬ 
tations  based  on  the  much  higher  speed  of  sound  of  the  water.  This  is  particularly  important  for 
plume  predictions  since  the  dynamics  of  interest  (formation  rise  and  fall  of  the  plume)  occur  on  the 
order  of  seconds.  The  disadvantage  of  an  incompressible  approach  is  that  the  early  time  shock 
effects,  in  particular,  the  early  time  dome  of  spray  ejected  upward  as  a  result  of  the  shock  reflecting 
off  the  surface,  is  not  predicted. 


COMPUTATIONAL  MODEL 

The  computational  model  used  in  both  the  two-  and  three-dimensional  codes  is  based  on  a  gen¬ 
eralized  formulation  of  hydrodynamics.9-12  This  formulation  uses  a  fixed  spatial  domain  Q,  where 
the  density  p,  velocity  u,  and  the  pressure  P  are  governed  by  the  mass  and  momentum  conservation 
equations 

P,  +  V-(pu)  =  0,  (2-la) 

(pu);  +  V-(puu)  =  -pg  k-VP ,  (2-lb) 

subject  to  the  constraint 

P^Po-  (2-lc) 

where  g  the  gravitational  constant,  -k  the  unit  vector  in  the  direction  of  the  gravitational  force,  and 

p0  is  the  constant  density  of  the  incompressible  liquid.  In  constructing  solutions  for  the  constrained 
system  (2-1),  it  is  assumed  that  liquid-on-liquid  collisions  behave  inelastically,  thereby  causing  a 
reduction  in  the  total  mechanical  energy  of  the  flow  field.  These  energy  losses  are  associated  with 
breakdowns  of  the  classical  theory  and  may  be  attributed  to  turbulence.13 

The  density  field  divides  Q,  into  two  time-varying  regions,  namely,  the  liquid  region  D(0  = 
{x  e  Q:  p(x,?)=Po}>  and  the  nonliquid  region  where  p  <  p0.  The  interfaces  separating  these  regions 
are  the  free  surfaces.  The  numerical  solution  “captures”  these  surfaces  as  slightly  smeared  inter¬ 
faces.  The  nonliquid  regions  are  characterized  by  specifying  uniform  pressures  in  each  of  its  con¬ 
nected  subsets.  For  example,  in  the  atmospheric  region  above  the  liquid  the  pressure  is  set  to  the 
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ambient  air  pressure.  In  an  underwater  bubble  the  pressure  PB  may  be  determined  using  the  adia¬ 
batic  law 

PB  =  CVJV,  (2-2) 

where  VB  is  the  instantaneous  bubble  volume,  C  is  a  constant,  and  y  is  the  ratio  of  specific  heats  of 
the  bubble  gases. 

Assume  that  the  density  and  velocity,  p"  ,u”  at  time  step  n ,  and  the  pressure  gradient  at  the 
previous  half  step,  VP”-1/2  are  known.  This  solution  is  evolved  using  the  following  three  step  time 
split  procedure. 

Convection 

The  solution  is  first  advanced 

(p\un)->(p,u) 

by  “solving”  the  conservation  laws  (2-la,b)  without  including  the  VP  term  on  the  right-hand  side 
of  (2- lb)  and  without  regard  to  the  constraint  (2-lc).  This  step  is  implemented  numerically  using  a 
formally  second-order  Godunov-type  method,  which  uses  slope  limiting  in  space  and  explicit 
predictor-corrector  time  stepping. 

Redistribution  of  Density  and  Momenta 

Next,  the  density  and  momenta  are  redistributed 

(P,u)->(p,u) 

so  that  the  constraint  (2-lc)  is  satisfied,  the  global  conservation  of  mass  and  momenta  are  main¬ 
tained,  and  the  energy  is  nonincreasing.  This  step  is  of  critical  importance  to  the  stability  and  con¬ 
vergence  of  the  overall  algorithm.14  The  density  is  redistributed  using  an  approximate  solution  to 
an  obstacle  problem.  This  solution  is  obtained  using  a  constrained  direction  preconditioned  conju¬ 
gate  gradient  method.  The  momenta  redistributions  are  determined  as  solutions  of  two  elliptic  self- 
adjoint  problems.  Discretizations  of  these  problems  yield  systems  of  linear  equations  with  diago¬ 
nally  dominant  matrices,  which  are  efficiently  solved  using  a  diagonally  preconditioned  conjugate 
gradient  method.  After  this  step,  p  =  pn+1  and  the  new  nonliquid  region  is  then  determined  along 
with  the  pressure  In  each  of  its  connected  subsets.  In  the  computational  space  the  new  liquid 
region,  D”+1  =  D(Jn+1),  is  defined  to  be  the  collection  of  grid  cells  Ctj  such  that 

P?jl  -  (l-eP)Po>  (2-3a) 

where  p”^1  is  the  density  in  cell  Ctj.  A  typical  value  used  for  ep  is  0.04. 

When  computing  shallow  depth  bubble  dynamics,  it  was  found  that  because  of  thin  layers  of 
water  between  the  bubble  and  the  air  surfaces  when  the  bubble  is  near  its  maximum  volume,  a 
value  £p  =  0.04  would  cause  the  bubble  to  “vent”  into  the  air  prematurely  unless  the  grid  size  was 
excessively  small.  This  situation  was  remedied  by  including  an  additional  factor,  zA .  All  cells  C,  ,, 
that  are  adjacent  to  “air”  cells  are  also  included  in  the  liquid  region  D”+1  if 
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Ptf  a  d-e^Po-  (2-3b) 

A  typical  value  used  for  zA  is  0.5.  Generally,  this  value  was  used  for  first  80  percent  of  the  dura¬ 
tion  of  the  first  bubble  cycle.  After  this  time,  the  value  of  zA  would  be  reduced  to  0.1.  This  was 
found  to  improve  stability,  as  monitored  by  computing  the  total  energy  at  each  time  step,  for  the 
long-time  bubble  dynamics. 

After  the  redistribution  step  u  =  un+1  in  the  new  nonliquid  region.  However,  u  is  not  con¬ 
sistent  with  (2- lb)  in  the  new  liquid  region. 


Pressure  Projection 

In  this  step  the  velocity  is  corrected 

u  u  rt+1 


using  the  gradient  of  the  new  pressure,  pn+l/l.  The  pressure  P  =  pn+Vl  is  the  solution  of 

TV- i-VP  =  V-u  in  D”+1  (2-4) 

pn+i  ' 

where  %  is  the  time  step.  This  equation  is  discretized  using  a  finite  element  method  with  bilinear 
elements  and  the  resulting  linear  system  is  solved  using  an  incomplete  Cholesky  preconditioned  con¬ 
jugate  gradient  method.  The  new  velocity 


n+i  _  xVP 
un  i  =  u  - 


n+Vi 


n+ 1 


is  divergence  free  in  Dn+1  and  thus  is  consistent  with  (2-lb). 


(2-5) 


To  determine  the  pressure  uniquely  using  (2-4),  boundary  conditions  must  be  specified.  On 
those  portions  of  the  boundary  of  D"+1  that  correspond  to  “wall”  boundaries  of  Q.  a  Neumann  con¬ 
dition  is  specified,  namely 


dP  n+l_ 
x—  =  p"+1«- 
an 


un. 


Note  that  this  condition,  together  with  (2-5),  implies  that  un+1n  =  0  along  wall  boundaries. 


Multiple  Bubble  Logic 

Dirichlet  conditions  on  the  pressure  are  specified  along  all  boundaries  of  the  nonliquid  regions. 
This  specification  is  complicated  by  the  fact  that  there  can  be  several  different  types  of  nonliquid 
regions  that  evolve  with  time.  As  mentioned  earlier,  the  nonliquid  regions  can  be  designated  as 
either  “air”  A  or  “bubble”  B  having  specified  pressures  PA  and  PB ,  respectively.  The  air  pressure 
PA  is  assumed  to  be  constant  for  all  time.  It  is  possible  that  the  bubble  region  may  become  discon¬ 
nected  and  in  general  we  have 

Kn 

B”  =  U** 

k= 1 
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where  Kn  is  the  number  of  distinct  bubbles  at  time  step  n ,  and  £"  denotes  the  disjoint  connected 
region  comprising  bubble  k  at  time  step  n .  Each  of  these  bubbles  can  have  a  different  pressure  j , 
and  as  long  as  these  regions  remain  distinct,  each  bubble  is  assumed  to  behave  adiabatically.  When 
a  bubble  splits  into  two  disjoint  regions,  say,  Bg  — »  £”+1  i^j  5"+1,  the  pressure  Pf+l  =  P”+1  is 
computed  assuming  the  volume  changes  occurred  before  the  split,  that  is, 


Pn+ 1  _  pre+1  _ 

l  ~Pm  ~ 


Similarly,  if  two  disjoint  bubbles  merge  into  a  single  connected  region,  say,  Bj1  {j  5”  — >  £/I+1  the 
new  pressure  is  given  by 

pn^  =  (?w + nvz,w  +  var1 

k  (y»+1)Y 

These  formulas  are  easily  extendible  to  general  cases  treating  any  finite  number  of  bubbles  merging 
and  splitting.11,12  Note  that  whenever  merging  occurs  the  pressure  of  the  new  connected  region 
changes  instantaneously.  Similarly,  whenever  a  bubble  comes  in  contact  with  the  air  region,  the 
bubble  is  said  to  “vent”  into  the  air  and  the  pressure  is  instantaneously  changed  into  the  air  pres¬ 
sure.  This  is  obviously  a  crude  approximation  to  what  is  physically  occurring  where  the  venting 
takes  a  finite  amount  of  time,  or  the  bubbles  can  partially  vent.  Nevertheless,  as  long  as  the  venting 
process  is  not  the  critical  phenomenon  being  studied,  we  can  still  hope  to  produce  reasonable  simu¬ 
lations  to  shallow  depth  explosion  bubbles. 

Another  capability  of  the  algorithm  is  the  ability  to  predict  the  inception  of  cavitation. 15  This  is 
accomplished  by  adding  another  constraint,  P  >  Pc  to  (2-1),  where  Pc  is  a  prescribed  cavitation 
pressure.  In  this  case,  (2-4)  is  replaced  by  an  obstacle  problem  for  the  pressure.  In  regions  where 
P  =  Pc ,  and  p  =  p0,  the  divergence  of  the  velocity  can  be  positive.  In  such  regions,  it  follows 
from  (2-la)  that  the  density  can  decrease.  When  this  occurs,  and  <  (l-ep)p0  the  cell  is  desig¬ 
nated  as  part  of  the  “cavitation”  region  C.  However,  because  of  truncation  error,  it  is  also  possible 
that  pi  j  <  (1— £p)p0,  while  P  >  Pc  in  the  interior  of  the  water  region.  This  is  referred  to  as 
“numerical  cavitation”  and  is  suppressed  by  simply  setting  p/'t1  =  (l-ep)p0  in  the  redistribution 
step.  Although  this  step  violates  conservation  of  mass,  it  reduces  the  truncation  error  in  the 
incompressible  water  region.  For  the  computations  presented  in  this  study,  the  total  percentage  of 
mass  added  per  bubble  period  was  only  about  0.005  percent. 
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CHAPTER  3 

BUBBLE  THEORY  AND  BENCHMARK  COMPUTATIONS 

This  chapter  presents  the  basic  theory  for  underwater  bubbles  based  on  an  incompressible 
liquid  assumption,  as  well  as  the  empirical  formulas  used  to  predict  bubbles  from  different  explo¬ 
sives.  The  empirical  formulas  are  examined  carefully  since  they  represent  approximations  that  are 
valid  for  moderately  deep  charges  only.  Some  background  information  on  the  properties  of  our  com¬ 
putational  algorithm  and  a  convergence  study  is  also  presented.  Theoretical  corrections  to  the  bubble 
period  in  the  proximity  of  the  free  surface  are  compared  to  computational  results  and  provide  an 
additional  validation  of  our  code  for  predicting  free  surface  effects.  Theory  is  presented  which 
shows  that  with  an  incompressible  water  model  and  an  adiabatic  gas  bubble,  the  maximum  volume 
of  the  bubble  will  be  reduced  as  the  bubble  is  moved  closer  to  the  free  surface  in  the  absence  of 
gravity.  An  empirical  formula  for  the  volume  reduction  is  provided  based  on  the  computational 
results.  The  theory  is  also  used  in  combination  with  bubble  period  and  maximum  size  measure¬ 
ments  to  determine  new  bubble  parameters.  These  new  parameters  are  more  consistent  with  shallow 
depth  observations,  than  parameters  derived  from  deeper  charges.  Finally,  a  comparison  to  an 
experiment  performed  by  Blake  and  Gibson4  is  used  as  a  validation  of  the  computational  predictions 
of  the  detailed  behavior  of  the  plume  formation  and  bubble  collapse. 


SPHERICAL  BUBBLE  THEORY 

The  equation  of  motion  for  a  spherical  bubble  in  an  infinite  incompressible  medium,  often 
referred  to  as  the  Rayleigh-Plesset  equation,1  is  given  by 

ad  +  -|a2  =  ~  PJ  (3-1) 

1  Po 


where  a  =  a (t)  is  the  bubble  radius  at  time  t,  PB  =  PB(t )  is  the  bubble  pressure  at  time  t,  is 
the  pressure  at  infinity,  and  p0  is  the  constant  density  of  the  incompressible  medium.  This  equation 
can  be  derived  from  the  conservation  equations  (2-la,b)  in  the  liquid  region.  Given  an  initial  radius 
A0  and  an  initial  bubble  pressure  PB,  the  pressure  at  later  times  can  be  derived  using  the  adiabatic 
assumption  (2-2), 


PB(t)  = 


Pb°(A°)3y 
(«(0)3y  ' 


Let 


a  = 


Lmax 
A  . 

^  mm 


(3-2) 
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be  the  ratio  of  maximum  to  minimum  bubble  radii.  If  initially  the  bubble  is  at  rest  (a  (0)  =  0)  and 
a  (0)  =  A  0  =  Amin,  an  integration  of  (3-1)  yields 


PB°  =  P  Jl-Y) 


1-a3 


l_a3(l-7) 

Solutions  to  (3-1)  are  periodic  and  the  period  Tx  can  be  determined  (e.g.,  Friedman16)  using 


t.  =  2/1(•^-r5'6(2Ilpor1,3 


Em 
p  5/6 


(3-3) 


(3-4a) 


where 


a  i 


W 


a 


3 /2da 


a  o 


V 1  -  a 3  -  ka  3^> 
A, 


L  =  A 


nun 


1  + 


a0  ~  ' 

p  o 

r B 


A  . 
^min 


PJy-1) 


1/3 


a  i  = 


*  = 


Lmax 


p  0 

*  ft 


P-iy-i) 


3y 


(3 -4b) 
(3-4c) 
(3-4d) 


and 


E  = 


4k  .  3 
2  ^min 


(3-4e) 


In  the  above,  a0  and  a  x  are  the  nondimensionalized  minimum  and  maximum  radii  and  correspond 
to  the  roots  of  the  denominator  in  the  integrand  of  (3-4b),  L  is  the  length  factor  and  E  is  the  total 
energy. 


NUMERICAL  ACCURACY  AND  CONVERGENCE 

The  first  set  of  benchmark  computations  are  comparisons  to  the  solution  of  a  spherically  sym¬ 
metric  bubble.  For  this  test,  the  code  BUB2D  was  used  with  the  acceleration  due  to  gravity  set  to 
zero,  together  with  the  assumption  that  the  bubble  gas  is  adiabatic  (y  =  1.3).  The  initial  conditions 
used  were 

P 

A0  =  1.0,  ~  =  342.864,  (3-5a) 

*  oo 

which  corresponds  to  the  “free  field”  value  Amax  =  10.  Furthermore,  we  used 

P<a=l0,  p0  =  0.031081,  (3-5b) 

so  that  the  free  field  period  is  =  1 .068407.  This  value  was  determined  by  both  using  an  adaptive 
Runge-Kutta  4-5  method17  with  a  tolerance  of  10-10,  on  problem  (3-1)  directly,  as  well  as  by  using 
(3-4a)  with  the  integral  (3-4b)  computed  using  the  midpoint  rule  with  repeated  Richardson 


3-2 


NSWCDD/TR-94/ 1 56 


extrapolation.  Since  the  integrand  has  a  singularity,  the  error  of  the  midpoint  rule  can  be  shown  to 

co 

be  of  the  form  Error  =  where  rt  =  i  -  Vi,  instead  of  ri  =  2 i  if  the  integrand  was  smooth. 

£=1 

This  problem  was  computed  in  axially  symmetric  (r,  z)  coordinates.  The  grids  consisted  of  a 
region  of  uniform  cells  of  size  Ar  =  Az  =  h  in  the  region  r  <  12,  and  Zc<z  <  Zc+12,  where  Zc  is 
the  z  location  of  the  bubble  center.  Symmetry  conditions  were  imposed  along  the  axis  r  =  0  and 
along  the  the  horizontal  line  through  the  center  of  the  bubble  Z  =  Zc .  The  grid  was  stretched  radi¬ 
ally  to  r  -  Rl,  and  upward  to  Zc  +  RL,  using  die  same  number  of  points  and  spacing  in  each  direc¬ 
tion.  The  time  steps  were  selected  using  an  adaptive  approach  based  on  changes  in  the  bubble 
volume.  The  number  of  steps  taken  until  the  first  period  T  (time  of  the  bubble  minimum  volume) 
is  denoted  by  NT.  The  value  ep  =  0.04  was  used  as  the  cutoff  parameter  (see  (2-3))  for  all  the  com¬ 
putations. 

A  summary  of  the  results  of  the  symmetric  bubble  simulations  are  displayed  in  Table  3-1. 

TABLE  3-1.  SYMMETRIC  BUBBLE  APPROXIMATIONS 


h 

Rl 

Nt 

Ah 

^max 

Ah 

"max 

Error  (%) 

j fh 

Error  (%) 

0.8 

50 

273 

9.312 

6.88 

0.9912 

7.23 

0.4 

100 

551 

9.608 

3.92 

1.0274 

3.84 

0.2 

200 

1100 

9.798 

2.02 

1.0479 

1.92 

0.1 

400 

2262 

9.902 

0.98 

1.0591 

0.87 

In  the  sequence  of  calculations  shown,  the  grid  was  refined  in  conjunction  with  increasing  the  loca¬ 
tion  of  the  boundary  RL,  as  well  as  increasing  the  number  of  time  steps  NT.  Notice  that  the  errors 
in  both  the  period  and  maximum  radius  are  approximately  halved  with  each  refinement  of  the  grid. 
This  is  indicative  of  first-order  convergence.  We  point  out  that  first-order  convergence  is  typical  of 
a  second-order  method  approximating  a  problem  with  a  discontinuity  (the  density  across  the  bubble 
surface).  A  first-order  method  would  typically  produce  a  convergence  rate  of  one-half  (i.e.,  two  lev¬ 
els  of  refinement  would  be  required  to  reduce  the  error  by  half).  These  results  can  also  be  used  to 
approximate  the  error  of  less  trivial  problems,  based  on  the  relative  grid  size  and  number  of  time 
steps  taken. 

EMPIRICAL  RELATIONS  AND  APPROXIMATIONS 

To  predict  underwater  explosion  behavior,  values  for  Amin,  a,  and  y  need  to  be  determined  as 
inputs.  Since  the  explosive  process  itself  is  not  being  modeled,  these  values  must  be  determined 
through  empirical  relationships. 

The  usual  assumption  for  the  energy  of  the  spherical  bubble  is  that  it  is  proportional  to  the 
charge  weight,  that  is 
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E  =  QW 


(3-6) 


where,  W  is  the  charge  weight  (mass)  and  Q  is  the  energy  per  unit  mass  released  by  the  explosive, 
after  the  shock  has  passed.16  It  follows  from  (3-6),  (3-4e),  (3-3),  and  (3-2)  that 


A  —  7 

^max  ~~  u  c 


l-cc3^ 


'll/s 


1-a  3y 


wm  wv3 

—  */  r 


1/3 


a  p  1/3 

*  an 


(3-7a) 


where 


= 


3  Q 


4% 


1/3 


(3-7b) 


is  constant  In  addition  to  (or  in  conjunction  with)  (3-6),  three  additional  empirical  relations 
describing  the  bubbles  produced  from  explosive  charges  are  often  cited: 1-3,18-25,29-32 

Amin  =  JNWm  =  qWll\  (3-8) 

1 1/3 

(3-9) 


and 


A/  ,max  ~  7 


tk  =k 


w 


w 


1/3 


.5/6 


(3-10) 


In  the  above,  J  is  the  charge  radius  “constant,”  N  is  the  charge  radius  ratio  “constant,”  K  is  the 
period  constant,  and  q  is  the  minimum  radius  constant.  Also,  P„  is  taken  to  be  the  hydrostatic 
pressure  at  the  charge  depth.  In  units  of  feet  and  feet  of  water, 


Px  =  d+PA  (3-11) 

where  d  is  the  initial  charge  depth,  and  PA  is  the  air  pressure.  For  sea  water 

Pa  ~  pAa  =  33,  feet  of  water .  (3-12) 

Equation  (3-8)  simply  reflects  the  assumption  that  the  initial  bubble  volume  is  proportional  to 
the  charge  weight.  Values  for  q  are  difficult  to  measure  photographically,  because  of  the  deforma¬ 
tion  of  the  bubble  upon  its  collapse.  This  deformation  is  reduced  for  very  deep  charges,  but  to  our 
knowledge,  it  has  been  directly  measured  for  spherical  Pentolite  charges23  (q  =  0.272±.005).  How¬ 
ever,  the  value  q  =  0.288  was  reported  by  Yennie  and  Arons22  for  trinitrotoluene  (TNT)  charges 
using  an  indirect  approach  based  on  extrapolating  a  curve  fit  to  bubble  radius-time  curves. 

Arons  et  al. 20-22  also  approximated  the  value  q  =  0.286  for  TNT  charges  at  a  depth  of  500  ft, 
using  pressure  pulse  values  and  incompressible  bubble  theory.  Snay  and  Christian24  improved  on 
this  procedure  by  using  the  difference  of  the  maximum  and  minimum  bubble  pulse  pressure  meas¬ 
ured  a  given  distance  from  the  bubble  center.  For  example,  if  AP  (2? )  =  Pmax(i?)  -  P^P)  is  the 
difference  in  the  maximum  to  minimum  values  of  the  pressure  at  a  distance  R  from  the  charge 
center  (occurring  at  the  time  of  the  first  bubble  minimum  and  maximum,  respectively)  then  it 
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follows  from  integrating  the  conservation  of  momentum  equation  together  with  (3-3)  (cf.  Reference 
24*)  that 


A P(R)  R 
P  A 

1  oo  ^max 


a3  -  1 
1  -  a3(1-7> 


-  a-1  +  1 


^  7) 


(3 -13a) 


and 


P  max(^  ) 


R 


=  (7-D 


Lmax 


aJ 


1 


1  -  a3(1^ 


(3-13b) 


Therefore,  given  measured  values  for  A P(R)  or  Pmax(i?)  -  Px,  and  Amax,  as  well  as  a  value  for  y, 
one  can  obtain  a  by  solving  (3-13a)  or  (3-13b)  for  a  and  determine  q  using  (3-2)  and  (3-8).  Com¬ 
putations  of  q  using  the  data  from  Arons  et  al. 20-22  are  listed  below  in  Table  3-2  for  TNT  charges 
using  both  y-  1 .25  and  y  =  1.3.  For  all  cases  the  pressure  was  measured  at  a  location 
R  -  Wmi 0.352,  and  the  maximum  radii  were  determined  using  the  approximate  formula  (3-9)  with 
J  =  12.6. 


TABLE  3-2.  DETERMINATION  OF  q  USING  PRESSURE  MEASUREMENTS  FOR  TNT 


W  (lbs) 

d  (ft) 

A P(R) 

Emax(R)  P co 

q  (Y=1.25) 

q  (y=1.30) 

Pco 

0.505 

250 

9.52 

8.88 

0.303 

0.321 

2.507 

250 

8.72 

8.09 

0.319 

0.338 

0.505 

500 

5.60 

5.10 

0.297 

0.314 

2.507 

500 

5.39 

4.89 

0.305 

0.322 

12.01 

500 

5.31 

4.80 

0.308 

0.325 

Average  Values 

0.306 

0.324 

These  values  were  obtained  by  inverting  W3  using  Newton’s  method  and  are  correct  to  the  number 
of  digits  displayed.  (We  actually  iterated  until  the  relative  difference  in  the  iterative  approximations 
were  less  than  10-14)  We  remark  that  the  same  values  for  q  were  obtained  using  (3-13b)  and  the 
data  for  Pmax(R )  -  P«,. 

Note  that  the  values  listed  are  higher  than  the  value  q  =  0.286  derived  by  Arons21  (using 
y  —  1.25)  and  reported  in  subsequent  publications.24’25  This  discrepancy  is  due  to  several  approxi¬ 
mations  performed  to  simplify  formula  (3- 13b)  in  Reference  21.  We  remark  that  in  addition  to  the 
assumption  that  the  bubble  gases  behave  adiabatically  throughout  the  bubble  oscillation,  the  addi¬ 
tional  assumption  that  the  bubble  remains  spherical  was  also  used  to  derive  (3-13).  It  is  this  latter 
assumption  which  is  most  suspect  as  the  phenomena  of  non-spherical  bubble  collapse  is  well  docu¬ 
mented.4,8-11  The  non-spherical  collapse  is  caused  by  gravity  effects  and  is  less  pronounced  for 

*  The  formulas  (3-13a)  and  (3-I3b)  are  equivalent  to  Equations  18  and  19  of  Reference  24.  This  can  be  seen  using 
Equations  11  and  14  from  that  same  reference. 
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smaller  charges  and/or  at  deeper  depths.  Furthermore  it  could  be  expected  that  a  non-spherical  col¬ 
lapse  would  cause  liquid  collisions  to  occur  (an  upward  moving  water  jet  would  collide  with  the 
upper  bubble  surface),  which  would  reduce  the  energy,14  causing  a  lower  pressure  at  the  minimum 
volume  than  would  occur  with  a  spherical  collapse.  This  phenomena  appears  to  be  indicated  in  the 
data  listed  in  Table  3-2  where  at  each  depth  the  smaller  charges  yielded  larger  values  for  the  pres¬ 
sure  differences. 

The  adiabatic  exponent  y  can  be  estimated  from  the  properties  of  the  explosion  by  products.25 
For  example,  if  the  percentages  by  weight  of  the  by-products  are  known,  then  the  adiabatic 
exponent  for  the  mixture  is  simply  the  ratio  of  the  weighted  sums  of  specific  heats.  From  the 
experimental  data  of  Jones  and  Miller26  as  well  as  Omellas27  the  by-products  of  TNT  yield 
y  =  1.29 —>  1.30.  This  is  slightly  higher  than  the  theoretical  result  of  Jones  and  Miller  who 
estimated  y  -  1.27.26  The  chemical  products  for  a  variety  of  explosive  types  have  been  computed 
by  Renner  and  Short28  using  an  empirical  model  for  the  equation  of  state  for  predicting  the  detona¬ 
tion  products.  Table  3-3  contains  a  list  of  charge  types  and  their  adiabatic  exponents  determined 
from  the  data  listed  in  Reference  28. 

TABLE  3-3.  VALUES  OF  y  FOR  VARIOUS  EXPLOSIVE  TYPES 


Explosive 

Density  ( glcm  3) 

7 

TNT  (cast) 

1.60 

1.29 

Pentolite 

1.67 

1.31 

COMP  C-4 

1.6 

1.34 

PBXN-103 

1.87 

1.23 

HBX-1 

1.72 

1.23 

In  equation  (3-9)  we  have  distinguished  the  approximation  for  the  maximum  radius  using  a 
charge  radius  constant  J  from  the  “exact”  formula  (3-7),  where  /  =  7a  is  not  a  constant.  These 
formulas  would  agree  if  the  ratio  a  was  a  constant  independent  of  the  depth.  However,  it  follows 
from  (3-2),  (3-7),  (3-8)  and  (3-1 1)  that  a  depends  on  the  depth  and  can  be  determined  using 

1 1/3 

=  Nx{d  +PA)m,  (3- 14a) 


1 


G(a),- 


l-a3^ 


l-of3Y 


where  N„  can  be  defined  using 


IV  = 


(3- 14b) 


In  practice,  if  q  is  known,  and  the  maximum  radius  is  measured  at  some  depth,  then  a  can  be  deter¬ 
mined  at  that  depth  using  (3-8)  and  (3-2).  Given  a  value  for  y,  a  value  for  J x  can  then  be  evaluated 
using  (3-7),  and  can  then  be  determined  using  (3-1 4b).  Figure  3-1  is  a  graph  of 
a  =  G~1(Nco(d+  33)1/3)  as  a  function  of  the  depth  d,  under  the  assumption  that  y  -  1.3  for  values  of 
N„  =  0.03,  0.025,  0.02  and  0.015.  These  represent  typical  values  for  Nx.  For  example,  using  the 


3-6 


N  SWCDD/TR-94/ 156 


data  presented  in  Reference  19  for  TNT,  the  average  value  of  Nm  =  0.0236  is  obtained  if  q  =  0.324, 
=  0.021  if  q  —  0.286.  Also  based  on  data  from  several  sources  for  Pentolite  charges  the  value 
=  0.0194  is  obtained.  From  this  figure  it  can  be  seen  that  for  a  typical  value  of  Nx  =  0.02,  a 
ranges  from  15.1  to  4.5  as  the  depth  varies  from  0  to  1000  ft.  It  can  easily  be  shown  that  a  ->  1  as 
d  — >  °o,  and  a  — » as  d  ->  -PA.  Perhaps  more  important  is  the  variation  of  Ja  =  Ja(d)  as  a 
function  of  d.  Figure  3-2  shows  the  ratio  J  JJa  as  a  function  of  depth  for  the  same  values  of  N^. 
At  shallow  depths  d  ~  0,  J„  is  approximately  3.1  percent  larger  than  Ja  for  charges  with 
Nm  =  0.02.  At  a  depth  of  300  ft  this  ratio  increases  to  6.8  percent  and  at  1000  ft  the  ratio  is  11.0 
percent.  This  means,  for  example,  that  a  value  for  J  determined  using  (3-9)  for  a  charge  at  a  depth 
of  300  ft,  would  underpredict  the  maximum  radius  of  a  shallow  charge  by  approximately  4  percent. 
A  measurement  for  /  based  on  a  charge  at  1000  ft  would  underpredict  a  shallow  charge  maximum 
radius  by  approximately  8  percent. 

Figure  3-3  displays  computed  values  for  both  J  using  (3-9)  and  using  (3-7)  for  charges  of 
TNT  at  various  depths.  The  data  used  for  this  figure  came  from  Reference  19,  Table  XV,  which 
contains  data  for  0.66-lb  charges  at  depths  between  300  and  600  ft,*  and  from  Table  XVIII,  which 
contains  data  for  56  and  295-lb  charges  at  depths  between  36  and  95  ft.  Notice  that  the  spread  in 
the  values  of  J  due  primarily  to  the  shallow  data.  The  average  values  were  J  =  12.6  for  the  deeper 
charges  and  J  =  13.1  for  the  shallow  shots.  However,  the  values  for  are  more  consistent,  having 
an  average  of  Jx  =  13.59  for  the  deeper  shots  and  /«,  =  13.68  for  the  shallow  depths  assuming 
q  =  0.286.  The  corresponding  deep  and  shallow  values  when  q  =  0.324  are  13.73,  and  13.75, 
respectively.  The  dependence  of  J  on  the  depth  has  been  known  and  reported  by  Snay,29  but  it  has 
been  ignored  since  for  the  cases  of  previous  interest  J  varied  by  about  4  percent  which  is  within  the 
experimental  accuracy. 

It  is  also  interesting  to  note  the  the  value  quoted  in  the  references  was  J  =  12.6  in  reports  prior 
to  i97852.20-22, 24-25, 29-31  while  J  =  13.1  was  listed  afterwards.3,32  This  time  dependence  on  the 
value  J  was  caused  by  simply  using  values  from  Reference  19  Table  XVTTT,  instead  of  Table  XV.33 
However,  since  values  for  J  for  other  charges  are  often  given  relative  to  values  of  TNT,  the  radius 
constants  for  all  types  of  charges  increased  by  1.04=  13.1/12.6  when  the  report  by  Swisdak32 
appeared.  These  problems  would  have  been  avoided  if  (3-7)  was  used  instead  of  (3-9). 

Increases  in  the  measured  value  of  J  were  also  reported  for  pressed  Tetryl  charges.19  At 
depths  greater  than  200  ft  the  average  value  of  J  =  12.8  was  reported  for  both  0.0558-lb  charges 
(Table  XVI)  and  1/2  lb  charges  (Table  XV).  However,  at  depths  between  1.75  and  5  ft,  the  average 
value  J  =  14.1  was  reported  in  Table  XVII. 

We  also  distinguish  the  period  TK,  computed  using  (3-10)  from  the  theoretical  period  Tx.  The 
bubble  period  can  be  measured  reliably  from  experiments  through  the  use  of  pressure  gauges. 
Measurements  for  the  period  have  been  plotted  as  a  function  of  the  hydrostatic  pressure  Px  for 

*  The  data  for  shot  G21F  was  excluded  since  Swift  and  Decius19  excluded  it  for  computing  the  average  value  of  the 
period  constant  K  (see  page  4  of  that  reference). 
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several  charge  types  and  have  been  shown  to  lie  very  near  the  empirical  relation  (3-10)  for  depths 
ranging  from  100  to  800  ft.1  However,  from  (3-4a)  and  assumption  (3-5),  the  formulas  (3-4a)  and 
(3-10)  are  not  equivalent  since  the  integral  Ix  depends  on  the  depth.  It  is  interesting  to  compare 
values  of  to  Tk-  We  can  do  this  by  comparing  measured  values  of  K  to  computed  values  of 

p  5/6 

Ka(d)  =  Ta 


W 


1/3 


in  Figure  3-4.  Both  sets  of  values,  (y  =  1.3,  q  =0.324,  =  13.73)  and  ( y=  1.3,  q  =0.286, 

Jx  =  13.60)  for  TNT  were  used  to  compute  using  (3-4).  Note  that  for  a  fixed  charge  weight 
W ,  and  a  specified  depth  d,  values  for  y,  q,  and  J ^  will  completely  determine  the  solution  to  (3-1). 
Therefore,  comparing  Ka(d )  to  measured  values  of  K  provides  a  consistency  check  for  the  theory. 
In  both  cases  the  theoretical  values  Ka(d)  are  slightly  larger  than  the  average  measured  value 
K  =  4.36  and  have  only  a  small  dependence  on  the  depth,  thereby  justifying  the  approximation  (3- 
10).  The  choice  q  =  0.286  gives  a  better  agreement  than  q  =  0.324.  The  theoretical  values  for 
Ka(d)  when  q  =  0.286  are  less  than  2  percent  larger  than  the  measured  values  at  depths  under 
800  ft.  This  is  in  reasonably  good  agreement,  especially  when  uncertainties  in  determinations  of  y, 
q ,  and  are  considered. 

Since  Ka(d )  is  completely  dependent  on  values  of  y ,  q,  and  that  is 


Ka(d)  =  n(y,  d), 


it  may  be  possible  to  approximate  one  of  these  constants  by  solving  II  =  K  for  the  parameter  in 
question.  Of  course,  II  is  a  non-linear  function  and  a  unique  solution  cannot  in  general  be 
guaranteed.  For  example,  suppose  that  y  =  y0  is  known  and  the  maximum  radius  Amax  has  been 
measured  for  a  charge  of  weight  W  =  W0  at  a  depth  of  d  =  dQ.  Furthermore  assume  that  the 
period,  and  hence  K  was  measured,  but  that  values  for  the  minimum  radius  and  the  pressure 
differences  are  not  available,  so  that  q  remains  unknown.  Since  by  (3-7)  the  value  for  /«,  depends 
on  y  and  a,  which,  in  turn,  depends  on  q ,  values  for  both  q  and  can  be  approximated  if  the 
equation 


n(y0,  q,J eo(7o.  q\  do)  —  k 

has  a  unique  solution  q  =  q0.  Figure  3-5  shows  the  dependence  of  II  as  a  function  of  q  for  the 
values  y0  =  1.2,  1.25,  1.3,  and  1.35.  Values  for  Amax  and  W0  were  taken  to  be  the  average  values 
for  TNT  reported  in  Reference  19  at  a  depth  of  approximately  dQ  =  300  ft.  At  this  depth  the  meas¬ 
ured  values  for  K  were  between  4.32  and  4.41  (shown  as  dashed  lines  in  Figure  3-5),  with  an  aver¬ 
age  of  approximately  4.36.  The  corresponding  values  of  q  lie  between  0.233  and  0.305,  with  the 
average  corresponding  to  0.265.  These  values  are  substantially  lower  than  the  value  q  =  0.324 
determined  using  (3-13)  based  on  the  data  in  Table  3-2.  As  mentioned  earlier,  this  discrepancy  may 
be  due  to  the  assumption  of  a  spherical  bubble  collapse,  which  was  used  to  derive  (3-13).  How¬ 
ever,  the  choice  q  =  0.286  yields  results  within  the  experimentally  measured  data.  The  correspond¬ 
ing  values  for  Jx  are  shown  in  Figure  3-6. 
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Values  for  J  can  been  determined  fairly  accurately  through  the  use  of  high  speed 
photography.18  Unfortunately,  for  many  of  the  newer  charge  types,  values  for  J  have  not  been 
determined  directly  from  photographs  due  to  the  extra  cost,  but  instead  were  inferred  using  a  value 
of  K  from  period  measurements.  It  has  been  reported30  that  values  of  J  IK  are  relatively  constant, 
ranging  from  J  IK  =  2.84  for  Pentolite,  to  J  IK  =  2.96  for  pressed  Tetryl.19  Therefore,  given  K,  one 
can  use  J  =  2.90 K  as  an  approximation,  if  photographic  measurements  are  not  available.  However, 
this  can  be  another  source  of  error  since  according  to  (3-7),  neither  J  nor  K  are  independent  of 
depth.  (Although  Figure  3-4  indicates  that  K  is  less  sensitive  than  J,  particularly  for  shallow 
depths).  However,  in  the  absence  of  more  reliable  data  the  approximation  =  3.1  K  can  be  used 
based  on  the  values  displayed  in  Figures  3-2,  3-5  and  3-6  for  TNT  discussed  above.  A  more  accu¬ 
rate  procedure  for  determining  would  be  to  solve  the  equation 

U(y,q,Jco,d)  =  K,  (3-15) 

assuming  that  y  and  q  are  known,  and  K  has  been  measured  at  some  depth  d. 

The  values  often  given  in  the  literature  for  the  radius  constant  J  were  determined  for  sea 
water,  while  the  tests  discussed  in  this  paper  are  for  fresh  water.  In  fresh  water  it  can  be  expected 
that  the  maximum  bubble  radius  will  be  larger  due  to  the  fact  that  fresh  water  is  less  dense  than  sea 
water.  The  correction  for  fresh  water  is30 

4,-  =  (3-16) 

where  the  ratio  33.9/33  is  the  ratio  of  air  pressure  in  units  of  feet  of  fresh  water  to  sea  water.  It  is 
also  the  approximate  density  ratio  of  sea  to  fresh  water.  In  the  case  of  fresh  water,  the  value 
pa  =pfa  =  33.9  should  be  used  with  JF  oa  in  (3-7).  (Actually,  the  value  34  was  used  instead  of 
33.9  in  Reference  30.) 

Two-Dimensional  Bubbles 

Suppose  several  identical  charges  are  placed  in  a  line  and  are  separated  by  a  distance  that  is 
considerably  less  than  their  free  field  maximum  radius.  Then  it  could  be  expected  that  this  “line  of 
charges”  could  be  approximated  (at  least  at  the  central  charge)  as  a  two-dimensional  problem  in 
which  the  bubble  is  cylindrical  and  has  infinite  length.  In  such  cases,  the  three-dimensional  problem 
with  multiple  charges  could  be  approximated  by  a  much  simpler  (and  computationally  less  inten¬ 
sive)  two-dimensional  problem. 

Let  each  individual  charge  have  mass  W,  and  be  separated  by  a  distance  S.  From  (3-8)  each 
of  these  charges  would  produce  a  spherical  bubble  having  a  minimum  (initial)  radius  given  by 

ASS)  = 

Let  A  ^ ’  denote  the  minimum  radius  of  the  cylindrical  bubble  used  in  the  two-dimensional 
approximation.  Imposing  the  condition  that  the  bubble  volume  per  unit  length  is  the  same  for  both 
the  multiple  spheres  and  the  single  cylinder  yields 
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|^e»)s=«(Ag)fs. 

Combining  this  with  (3-8)  yields 

A^)=7=(A<g,)3'2  =  ^|«3'2M,n  =  ?(2D)M1'2,  (3-1 7a) 

where 


4<2D>  =  -|-93'2, 


(3- 17b) 


and  M  ~W!S  is  the  mass  per  unit  length  of  the  line  charge.  The  initial  pressure  for  the  cylindrical 
bubble  is  the  same  as  for  each  individual  charge  given  in  (3-3)  with  a  =  a(3D)  being  the  ratio  of  the 
individual  spherical  bubble’s  maximum  to  minimum  radii.  In  practice,  the  value  for  oS3D)  could  be 
determined  by  solving  (3- 14a).  An  estimate  for  the  maximum  radius  of  the  cylindrical  bubble  can 
be  obtained  by  equating  the  equivalent  expression  for  the  initial  pressure 


=  ^(1-7) 


l-(a(2P))2 

(l-(a(20))2M 


where  a(2D)  =  to  (3-3)  which  yields 


a(2D)  =  (a(3Z»)3/2 


Using  this  with  (3-17)  and  (3-7)  we  obtain 


A  (2D)  _  j  (2D) 
^  max  °  oo 


l_(a(2Z?)r2y  J 


1/2 


M 


1/2 


1/2 


(3- 18a) 


where 

j£D)  =  (3-1 8b) 

Therefore,  the  empirical  constants  and  initial  conditions  for  a  cylindrical  bubble,  can  be  determined 
directly  from  the  spherically  symmetric  bubble  theory. 


FREE  SURFACE  EFFECTS 

Here  the  effects  of  the  free  surface  on  the  period  and  maximum  bubble  volume  are  examined. 
Both  theoretical  and  computational  results  are  presented.  This  analysis  is  crucial  for  a  better  under¬ 
standing  of  the  comparisons  of  our  computational  results  to  the  experimental  data  presented  in  the 
next  chapter. 

Free  Surface  Effects  on  Bubble  Period 

The  effects  of  both  the  surface  and  a  rigid  bottom  boundary  on  the  bubble  period  has  been  stu¬ 
died  by  Friedman.16  Under  the  assumption  that  the  bubble  remains  spherical,  as  well  as  the  usual 
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assumptions  that  the  bubble  gases  behave  adiabatically  (2-2)  and  are  surrounded  by  a  constant  den¬ 
sity  incompressible  liquid,  he  used  potential  theory  together  with  the  “method  of  images”  to  derive 
a  formula  for  surface  and  boundary  corrections.  This  formula  can  be  expressed  as 


kT^PJ?{x) 

d 


(3-19) 


where  b  is  the  initial  distance  from  the  bubble  center  to  the  bottom  surface,  Td  b  is  the  corrected 
period,  is  the  “free  field”  period,  F(x)  is  defined  by  the  infinite  series 


where 


F{ *)  = 


x 

1  +  X 


+  (l-x) 


f.  c-iy* 

i  (2n+l)2  -  x2 


(3-20a) 


*  =  JTl’  ^ 

and  k  is  a  constant  depending  on  the  adiabatic  exponent  y  and  the  minimum  to  maximum  ratio  a 
(3-2).  In  particular,  if  values  of  Amin  and  Amax  are  known  a  priori  then  k  can  be  computed  using 


k  =  1.158 — -P 

/? 


(3-21) 


where  I{  is  defined  in  (3-4b)  and 


_ a5,2da _ 


(3-22) 


In  (3-21)  the  coefficient  1.158  =  0.25V2g/3,  where  g  =  32.174  is  the  acceleration  due  to  gravity  in 
units  of  ftlsec2.  Note  that  if  the  water  is  infinitely  deep,  b  =  °o,  it  follows  from  (3-20)  that  x  =  1 
and  F(x  )  =  Vi.  For  details  on  the  derivations  of  these  formulae  see  Friedman16  and  Cole.1* 

Figure  3-7  displays  the  comparison  between  the  theoretical  correction  TdlT^  (using  (3-19))  and 
the  computations  using  the  medium  grid  (h=0,2).  The  computational  results  plotted  are  Tj/T* , 
where  it  =  1.0479  (from  Table  3-1)  was  used  instead  of  =  1.0684  to  reduce  the  effects  of 
discretization  errors.  Using  (3-22)  together  with  the  initial  conditions  (3-5)  yields  a  value  of 
k  =  1.189,  which  is  typical  of  values  which  have  been  reported  (1.13  <  k  <  1.22)  from  experimental 
data.1*16’25  We  remark  that  formula  (3-19)  has  been  used  in  practice  in  its  approximate  form,  with 
Tk  replacing  T„  on  the  right  hand  side.25  Let  C  to  be  the  scaled  depth, 


A 

■^max 


(3-23) 


*  The  formula  (3-20)  is  taken  directly  from  Friedman,16  page  190.  Cole1  uses  a  different  equivalent  formula  but  er¬ 
roneously  omits  the  alternating  sign  (-1)”  whenever  dealing  with  formulas  resulting  from  the  image  theory  (page  323, 
330,  and  formula  (8.60)  on  page  333).  However,  his  Figure  8.1  on  page  335  for  die  image  function  is  correct.  That  is, 
it  represents  values  not  from  his  formula  (8.60)  but  values  which  are  obtained  using  an  alternating  series  in  (8.60). 
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Figure  3-7  shows  the  agreement  between  the  computations  and  theory  to  be  very  good,  particularly 
for  values  of  C>  1.5.  However,  the  disagreement  increases  for  smaller  values  of  C.  In  this  case, 
the  validity  of  the  theory  becomes  questionable  due  to  the  assumption  that  the  bubble  remains 
spherical  throughout  its  entire  motion.  An  additional  assumption  used  for  the  derivation  of  (3-19)  is 
that  the  boundary  condition  for  the  potential  at  the  free  surface  is  linear  and  homogeneous.  It  has 
been  more  recently  demonstrated34  that  the  more  correct  nonlinear  terms  in  the  free  surface  boun¬ 
dary  condition  become  extremely  important  for  C  <  1.5.  The  actual  computational  results  are  also 
listed  below  in  Table  3-4  together  with  the  computed  equivalent  maximum  bubble  radii. 

TABLE  3-4.  FREE  SURFACE  EFFECTS  ON  PERIOD  AND  MAXIMUM  RADIUS 


c  d 

Thd 

•4  max  (U  ) 

C~  A 

**  max 

0.56 

0.7655 

9.384 

0.60 

0.7759 

9.460 

0.80 

0.8200 

9.666 

1.00 

0.8554 

9.736 

2.00 

0.9474 

9.789 

4.00 

1.0020 

9.796 

oo 

1.0479 

9.797 

For  these  computations,  a  uniform  grid  with  h- 0.2  was  used  in  the  region  0  <  r  <  12,  and 
~{d+ 12)  <  z <12,  with  grid  stretching  to  r  -  RL  =  200,  and  downward  to  z  =  ~(d+RL).  The  value 
eA  =  0.4  was  used  at  the  air-water  interface  and  ep  =  0.04  at  the  bubble  interface  (see  Chapter  2). 

Free  Surface  Effects  on  Maximum  Radius 

Very  little  attention  has  been  given  to  the  effects  on  the  maximum  bubble  radius  due  to  the 
proximity  of  the  surface.  Here,  by  bubble  “radius”  we  mean  the  equivalent  radius  of  a  sphere  with 
the  same  volume,  since  in  general,  the  bubble  may  diverge  from  its  spherical  shape  at  the  time  of  its 
maximum  volume.  The  formulas  for  maximum  radius  found  in  Friedman16  are  not  appropriate 
since  their  derivation  depended  on  the  additional  assumption  that  the  velocity  of  the  surface  is  negli¬ 
gible  (Ibis  is  tbe  assumption  ^hich  linearizes  the  free  surface  boundary  condition)  and  the  bubble 
stays  spherical.  Indeed,  those  formulas  for  maximum  radius  predict  a  reduction  due  only  to  the 
vertical  migration  of  the  bubble  and  are  not  appropriate  here. 

Let  D(f )  be  a  domain  occupied  by  an  incompressible  liquid,  that  is, 

D(0  =  jx  :  p(x,t)  =  p0,  V-u(x,f)  =  oj. 

Consider  the  following  class  of  free  surface  problems. 
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DEFINITION:  ( Class  A  Problems ) 

A  Class  A  free  surface  problem  is  one  which  satisfies  the  following  assumptions: 

(i)  If  D(f)  extends  to  infinity  in  one  or  more  directions  then  the  velocity  and  pressure  have 
the  asymptotic  behavior 

u(x,/)  =  O 
P(x,t)  =  Poa  +  0 

as  1x1  — »  co,  where  P„  is  a  constant. 

(ii)  The  liquid  region,  D(t)  completely  surrounds  a  bubble  region  B(t ),  having  a  boundary 
denoted  by  dB(t).  This  region  contains  a  uniform  pressure  given  by 

PB=CVB-l  =  f(yB),  (3-24) 

where  C  >  0,  and  7  >  1  are  constants  and  VB  is  the  volume  of  B  (?)  given  by 

v„(f )  =  r  dx. 

B(t) 

Furthermore,  the  pressure  is  assumed  to  be  continuous  so  that  PB  =f(VB)  on  dB(t)  also. 

(iii)  All  other  finite  boundaries  of  D  are  such  that  either  u-n  =  0  (rigid  wall  or  symmetry 
boundaries)  or  P  =  Px  (free  surface). 

(iv)  No  gravity  is  acting. 

(v)  Initially,  the  liquid  has  zero  velocity,  that  is, 

u(x,0)  =  0.  (3-25) 

The  simplest  example  of  a  Class  A  problem  is  the  spherically  symmetric  bubble  in  an  infinite 
domain.  This  problem,  governed  by  equation  (3-1),  yields  periodic  solutions  for  the  bubble  radius 
with  Am}n  <a(t)<  Amax. 

Let 

^min  —  ^max  =  max  >  (3-26) 

denote  the  minimum  and  maximum  volumes  of  a  spherically  symmetric  bubble. 

Without  loss  of  generality  we  can  assume  that  initial  conditions  for  our  Class  A  problems 
satisfy 

=  (3-27) 

and 

PB(0)  =  Pg>Px.  (3_28) 
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Let  e{t)  be  the  kinetic  energy  of  the  liquid 

e(t)  =  [  V^plu \2dx. 
D(t) 


(3-29) 


The  following  result  states  that  an  upper  bound  on  the  bubble  volume  is  the  maximum  free  field 
bubble  volume. 

THEOREM  3-1:  The  volume  of  a  bubble  B(t)  governed  by  a  Class  A  free  surface  problem,  with 
initial  conditions  (3-27),  (3-28)  satisfies  the  inequality 

^min  ^  VB(t )  <  Vsmax,  for  t  >  0,  (3-30) 

with  equality  holding  only  if  e{t)  =  0. 

PROOF:  The  energy  balance  equation  is 

+  J  (JP  ~  PJ(*n)dS  =  0,  (3-31) 

dt  35(f) 

where  n  is  the  unit  outer  normal  to  D(t),  and  e  is  defined  by  (3-29).  This  is  derived  by  taking  the 
dot  product  of  the  velocity  vector  with  the  momentum  equations  (2- lb)  (with  g=0  by  Class  A  con¬ 
dition  (iv)),  integrating  in  space,  and  applying  the  Class  A  assumptions  (i)  and  (iii)  together  with  the 
Divergence  theorem.  Since  dB  (t)  is  the  surface  of  the  bubble,  its  speed  of  propagation  in  the  direc¬ 
tion  of  n  is  given  by  -u-n,  it  follows  that 

dV g 

(u-n)  =  — (3-32) 


Using  (3-32)  in  (3-31)  it  follows  that 

de 

dt 

or  equivalently 


UVfi 

f(yB)-p„-jr  =  o, 


e  +  F  (VB)  —  c , 


where  C  is  a  constant  of  integration,  and 


FOO  =  -J[f(t)-PJdT 


Note  that  for  V  >0 


and  using  (3-24), 


dF 

dV 


=  -[f(V)-P  J, 


S  =  ^(V)>0’ 

so  that  F  is  a  convex  function  of  V,  with  a  unique  minimum 


(3-33) 

(3-34) 

(3-35) 
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F(VJ  =  min  {F(V)>, 
v>o 

where  /(FJ  =  by  (3-35).  The  Class  A  assumption  (v)  implies  that  e(0)  =  0.  This,  combined 
with  (3-33)  and  the  initial  condition  (3-27)  yields 

e(t)=F(V*Din)-F(yB(t)). 

Noting  that  F(F^iax)  =  F (V ^in )  (from  (3-3)),  and  e(t)  >  0,  the  theorem  follows  from  the  convexity 
of  F. 


Theorem  3-1  states  that  the  presence  of  the  free  surface,  or  a  rigid  boundary,  will  in  general 
reduce  the  maximum  bubble  volume.  The  same  conclusion  holds  if  the  initial  shape  of  the  bubble 
is  non  spherical.  A  simple  explanation,  which  summarizes  the  proof  of  the  theorem,  is  that  if  there 
is  any  kinetic  energy  in  the  incompressible  liquid  at  the  time  of  the  maximum  bubble  volume  is 
attained,  then  the  potential  energy  of  the  bubble  F,  and  hence  volume,  must  be  less  than  if  the 
kinetic  energy  were  zero,  as  it  is  in  the  spherically  symmetric  case. 


Figure  3-8  is  a  graph  of  the  reduction  of  the  equivalent  maximum  radius  Amax(C)/Amax(°o), 


using  both  the  computational  results  from  Table  3-4,  as  well  as  an  empirical  fit  to  this  data,  given 
by  the  formula 


Amax(^~ ) 

A  max(°°) 


,  0.0071 

r  o  J 

L. 

0.007 1 A  ^(qq) 
d 3 


(3-36) 


If  a  measured  value  for  the  maximum  radius  is  used  in  Amax(C),  then  Amax(°o)  can  be  approximated 
by  solving  (3-36).  Of  course,  this  approximation  may  only  be  valid  when  a~10,  but  this  is  a  fairly 
typical  value  according  to  Figure  3-1. 


NEW  BUBBLE  PARAMETERS 

As  was  demonstrated  earlier  in  this  chapter,  the  commonly  used  empirical  formula  (3-9)  is 
inaccurate  for  shallow  charge  depths.  Therefore,  the  value  defined  in  (3-7)  needs  to  be  deter¬ 
mined,  as  opposed  to  J =  Ja,  which  is  depth  dependent.  First,  consider  direct  observations  of  the 
maximum  bubble  radius  for  shallow  charges  of  C-4,  and  both  shallow  and  deep  charges  of  Pentolite. 
Table  3-5  lists  the  computed  values  for  /«,  based  on  measured  maximum  bubble  radii  for  C-4  and 
Pentolite  charges.  In  this  table,  the  data  sources  prefixed  with  a  CD  describes  measurements  taken 
from  the  NSWCCD  Test  Pond  in  May,  1993. 
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TABLE  3-5.  MAXIMUM  RADIUS  MEASUREMENTS  AND  J„  VALUES 
FOR  PENTOLITE  AND  C-4  CHARGES 


Data 

Source 

Charge 

Type 

d 

(ft) 

W 

(lb) 

A  meas 
Amax 

(ft) 

4  max(°°) 

(ft) 

C 

3  a 

CD-9 

C-4 

2.5 

1.25 

4.57 

4.78 

0.52 

14.7 

15.0 

CD-11 

C-4 

5.0 

1.25 

4.75 

4.78 

1.05 

15.0 

15.3 

CD-17 

C-4 

9.0 

1.06 

4.08 

4.08 

2.21 

14.0 

14.3 

CD-10 

Pentolite 

2.5 

1.72 

5.50 

5.95 

0.42 

16.5 

16.8 

CD-13 

Pentolite 

5.0 

1.73 

5.75 

5.81 

0.86 

16.4 

16.8 

CD-16 

Pentolite 

9.0 

1.73 

4.84 

4:85 

1.86 

14.1 

14.5 

Ref.  19 

Pentolite 

335 

0.551 

1.46 

1.46 

229 

12.7 

13.7 

Ref.  19 

Pentolite 

396 

0.548 

1.33 

1.33 

297 

12.2 

13.2 

Ref.  23 

Pentolite 

5982* 

1.06 

0.698 

0.698 

8570 

12.4 

15.0 

Ref.  23 

Pentolite 

11852f 

1.03 

0.545 

0.545 

21746 

12.3 

15.6 

Obviously,  there  is  a  fairly  large  amount  of  scatter  in  the  data  listed  above.  In  addition  to 
modeling  assumptions,  the  scatter  can  also  be  attributed  to  inaccuracies  in  measuring  the  bubble 
size,  and  variations  of  the  explosive  product  itself  (e.g.,  age,  storage  conditions,  and  manufacturing 
process  may  all  have  some  effects  on  the  detonation.)  However,  it  appears  that  the  values  for  Jm  at 
shallower  depths  are  somewhat  larger  than  the  values  determined  at  greater  depths  with  the  excep¬ 
tion  of  the  very  deep  shots  reported  in  Reference  23.  This  discrepancy  cannot  be  explained  by  our 
theory.  Indeed,  equation  (3-36)  and  Theorem  3-1  predict  the  opposite  effect.  Therefore,  if  this 
effect  is  real  then  it  must  be  due  to  other  physical  effects  (shock  effects,  dissolved  gases  near  the 
free  surface,  etc.)  not  modeled  in  our  theory.  The  result  for  the  d  =  2.5  Pentolite  case  (CD- 10)  may 
also  not  be  valid  here  since  it  is  likely  that  the  bubble  may  have  partially  vented  into  the  atmo¬ 
sphere  before  it  attained  its  maximum  size. 

As  mentioned  previously,  an  alternative  method  for  determining  J  x  is  to  use  period  measure¬ 
ments  and  solve  equation  (3-15).  Since  the  free  surface  and  bottom  locations  had  an  effect  on  many 
of  the  tests  reported,  the  values  J  x  were  determined  using  (3-15)  with 

p  5/6 

n(y,  q ,  Joo,  d)  —  “yj 

with  Tb  d  defined  using  (3-19).  The  tabulated  results  of  this  procedure  are  listed  in  Table  3-6 
below.  In  this  table,  the  values  listed  for  Amax(°o)  were  determined  from  the  derived  values  of 
from  solving  (3-15)  and  using  (3-7a),  with  a  determined  from  (3-14).  Here,  data  sources  prefixed 
with  ARV  describe  tests  performed  in  Arvonia  in  August,  1993,  and  the  BP  prefix  corresponds  to 
tests  at  the  Briar  Point  facility  in  June,  1994. 

*  Average  values  for  3  shots  at  depth  between  5934  and  6030  ft. 
t  Average  values  for  9  shots  at  depths  between  11812  and  11885  ft. 
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TABLE  3-6.  COMPUTED  Jm  VALUES  BASED  ON  PERIOD  MEASUREMENTS 


Data 

Source 

Charge 

Type 

w  ■ 

(lbs.) 

d 

(ft.) 

b+d 

(ft.) 

j'meas 

(sec. ) 

^max(°°) 

(ft.) 

C 

CD-11 

C-4 

1.25 

5.0 

24 

0.1952 

4.48 

1.12 

14.4 

CD-3 

C-4 

1.25 

9.0 

24 

0.1952 

4.21 

2.14 

14.0 

CD-17 

C-4 

1.06 

9.0 

24 

0.1861 

3.99 

2.25 

14.0 

ARV-1 

C-4 

10.0 

6.0 

110 

0.34 

9.63 

0.62 

15.6 

CD- 13 

Pentolite 

1.75 

5.0 

24 

0.2120 

4.95 

1.01 

14.3 

CD-4 

Pentolite 

1.75 

9.0 

24 

0.2150 

4.67 

1.93 

14.0 

CD- 16 

Pentolite 

1.72 

9.0 

24 

0.2157 

4.68 

1.92 

14.1 

Ref.  35 

Pentolite 

0.671 

39 

87 

0.108 

2.78 

14.0 

13.8 

Ref.  35 

Pentolite 

0.655 

39 

106 

0.108 

2.79 

14.0 

14.0 

Ref.  19 

Pentolite 

0.552 

335 

3600 

0.02607 

1.46 

229 

13.6 

Ref.  19 

Pentolite 

0.566 

404 

4800 

0.02290 

1.37 

295 

13.6 

Ref.  19 

Pentolite 

0.550 

396 

4700 

0.02240 

1.32 

300 

13.2 

Ref.  19 

Pentolite 

0.551 

380 

6700 

0.02347 

1.36 

279 

13.4 

BP-7 

HBX-1 

206 

25 

54 

0.883 

23.18 

1.08 

16.4 

BP-4 

HBX-1 

206 

26 

54 

0.894 

23.34 

1.11 

16.6 

BP-6 

HBX-1 

206 

27 

54 

0.891 

23.21 

1.16 

16.6 

BP-2 

HBX-1 

206 

28 

54 

0.891 

23.01 

1.22 

16.6 

BP-5 

HBX-1 

206 

29 

54 

0.889 

22.80 

1.27 

16.5 

BP-3 

HBX-1 

206 

30 

54 

0.888 

22.67 

1.32 

16.5 

BP-8 

PBXN-103 

26.2 

10.1 

54 

0.542 

15.08 

0.67 

19.1 

BP-9 

PBXN-103 

26.2 

13.8 

54 

0.544 

13.69 

1.01 

17.9 

Ref.  35 

PBXN-103 

29.3 

30 

65 

0.550 

13.88 

2.16 

19.3 

As  in  Table  3-5,  the  estimations  for  appear  to  be  larger  for  smaller  values  of  C  (with  the  excep¬ 
tion  of  the  PBXN-103  charge  at  d  =  30).  In  this  case,  however,  this  can  be  attributed  to  the 
underestimation  of  the  period  using  (3-19)  (see  Figure  3-7).  Also,  the  values  obtained  for  using 
the  period  are  substantially  below  the  corresponding  values  listed  in  Table  3-5  for  both  C-4  and 
Pentolite  charges  with  d  =  5  and  d  =  9.  The  values  obtained  for  Pentolite  charges  at  depths  of 
d  =  335  and  d  =  396  using  both  procedures  yield  essentially  the  same  result.  Unfortunately,  other 
comparisons  can  not  be  made  at  this  time  since  the  periods  for  the  Pentolite  charges  listed  in  Refer¬ 
ence  23  were  not  listed,  and  the  bubbles  were  not  photographed  for  the  HBX-1  or  PBXN-103  shots 
listed  in  Table  3-6. 

The  explosive  bubble  parameters  used  in  the  next  chapter  are  listed  below.  The  2D  values 
listed  were  derived  using  (3-17)  and  (3-18). 
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TABLE  3-7.  EMPIRICAL  BUBBLE  PARAMETERS 


Charge 

7 

9 

q(2D) 

C-4 

1.34 

0.286 

15.3 

0.177 

69.1 

Pentolite 

1.30 

0.272 

15.3 

0.164 

69.1 

HBX-1 

1.23 

0.360 

16.5 

0.249 

77.4 

PBXN-103 

1.23 

0.400 

18.8 

0.292 

94.1 

These  values  for  C-4  and  Pentolite  are  only  roughly  based  on  the  data  listed  in  Tables  3-5  and  3-6. 
For  C-4  the  value  J ^  =  15.3  is  simply  the  value  in  the  case  when  C  ~  1.0  (shot  CD-11)  based  on 
photographic  measurement.  It  also  corresponds  to  the  shallowest  case  for  which  the  bubble 
remained  nearly  spherical  at  the  time  of  its  maximum  size.  Due  to  the  wide  scatter  in  the  values 
listed  in  Table  3-5  for  Pentolite  we  simply  chose  the  same  value  as  C-4.  The  equivalence  of  C-4 
and  Pentolite  is  supported  by  tests  performed  by  Heathcote  and  Niffenegger.36  These  values  should 
be  considered  valid  for  shallow  depths  only,  and  Jx  ~  13.8  is  probably  a  better  approximation  for 
C  >2.  Note  that  this  value  is  much  closer  to  the  value  J ^  —  13.6  for  TNT  charges  derived  earlier. 
For  HBX-1  and  PBXN-103  the  averages  of  the  values  for  listed  in  Table  3-6  was  used. 

SPARK  BUBBLE  COMPARISONS 

We  continue  our  code  validations  by  attempting  to  match  the  excellent  photographs  of  a  shal¬ 
low  depth  spark  generated  bubbles.4  These  experiments  were  done  in  free-fall,  thereby  removing  the 
effects  of  gravity.  Furthermore,  it  can  be  expected  that  the  shock  wave  from  the  spark  will  be 
weaker  than  that  resulting  from  an  explosive  so  the  early  time  compressible  effects  (e.g.,  cavitation 
and  spalling  at  the  surface)  are  not  evidenced  in  these  experiments. 

In  particular,  consider  the  case  C  =  0.56,  In  the  experiments  the  value  Amax  =  Dmax/2  was 
used  as  a  value  for  the  maximum  radius  for  the  determination  of  the  scaled  depth  (3-23),  where 
Dmax  was  the  maximum  observed  horizontal  diameter  of  the  bubble.  According  to  the  computations 
we  conducted,  the  maximum  width  matched  the  maximum  equivalent  radius  to  less  than  0.5  percent 
for  examples  with  0.5  <  C  <  0.6  despite  the  fact  that  the  bubble  does  not  remain  spherical.  (By 
maximum  equivalent  radius  we  mean  the  radius  of  a  sphere  with  the  same  volume  as  the  maximum 
bubble  volume.) 

In  order  to  slightly  improve  our  simulation  of  the  experiment  we  used  the  free  surface  correc¬ 
tion  based  on  (3-36).  In  particular,  for  the  initial  conditions  we  modified  (3-5),  using  A0  =  1.042 
instead  of  A0  =  1.0,  so  that  A^C)  =  DmJ2  ~  10,  when  the  initial  depth  is  d  =  5.6.  The  grids 
used  for  the  computations  consisted  of  a  region  of  uniform  cells  of  size  h  in  the  region  0  <  r  <  12 
and  -16  <  z  <  24.  Outside  of  this  uniform  grid  region  the  grids  were  stretched  downward  to 
z  =  Zb  =-119,  upward  to  z  =  Zt  =  100,  and  radially  to  r  =  RL  =72.6,  which  corresponds  to  a 
scaled  version  of  the  experimental  apparatus. 

Figure  3-9  compares  a  copy  of  the  photographs  of  the  experiments  of  Blake  and  Gibson  to  the 
computations  using  the  “fine”  grid  ( h  =0.1,  160  x  488  cells,  Nr  ~  2000  steps).  The  times  printed 
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above  the  figures  are  scaled  using 


where  t  is  the  actual  time.  The  frames  corresponding  to  the  computational  results  contain  density 
contours  of  p  =  0.5p0  at  the  same  scaled  times  as  the  photographs.  Note  the  agreement  in  predict¬ 
ing  the  large  central  jet  rising  above  the  surface,  as  well  as  the  smaller  downward  moving  jet,  which 
can  be  seen  in  the  frames  corresponding  to  T  =  0.59.  This  downward  moving  jet  passes  through 
the  center  of  the  bubble  and  by  T  =  1.16  has  impacted  its  bottom  surface.  This  jet  impacts  the  bot¬ 
tom  slightly  earlier  in  the  computations,  but  this  may  be  due  to  the  fact  that  the  gases  in  the  spark 
bubble  are  very  likely  at  or  near  the  vapor  pressure  Pv  as  opposed  to  having  an  adiabatic  behavior. 
In  this  event  the  times  should  be  scaled  using  PM  -  Pv  in  place  of  Pm  in  (3-37). 

Figure  3-10  shows  the  results  of  computations  using  both  the  “medium”  (h  =  0.2,  80  x  264 
cells,  Nt  ~  1000  steps)  and  “coarse”  ( h  =  0.4,  41  x  141  cells,  NT  ~  500  steps)  grids.  This  shows 
that  qualitative  agreement  can  be  attained  without  large  computational  expense.  Quantitatively,  the 
scaled  heights  of  the  plume  from  both  the  experiments  and  the  computations  are  shown  in 
Figure  3-11.  The  computed  plume  heights  are  all  above  the  measured  values,  which  again  may  be 
due  to  the  vapor  pressure.  Nevertheless,  the  agreement  is  still  very  good.  The  run  times  for  these 
computations  were  approximately  15  min  for  the  coarse  grid  run,  3  hr  for  the  medium  grid  run,  and 
24  hr  for  the  fine  grid  run  on  an  HP  9000-735  workstation. 
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FIGURE  3-1.  DEPENDENCE  OF  THE  RADIUS  RATIO  ON  THE  DEPTH 
WITH  y=  1.3  AND  VARIOUS  VALUES  OF  N„ 


d  (ft.) 

FIGURE  3-2.  RATIO  BETWEEN  THE  RADIUS  CONSTANT  J x  AND  Ja(d) 
WITH  7  =  1.3  FOR  VARIOUS  VALUES  OF  N„ 


3-20 


J  and  Values 


NSWCDD/TR-94/156 


1.35 

1.3 

1.25 

1.20 


K  (Measured)  - 


FIGURE  3-5.  THEORETICALLY  DETERMINED  VALUES  OF  Ka  AS 
A  FUNCTION  OF  q  FOR  VARIOUS  VALUES  OF  y 


FIGURE  3-6.  DEPENDENCE  OF  ON  q  FOR  VARIOUS  VALUES  OF  y 
(BASED  ON  A  VALUE  Amax  MEASURED  AT  d  =  300  FT) 
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FIGURE  3-7.  COMPARISON  OF  THEORY  TO  COMPUTATIONS  OF  FREE  SURFACE 
EFFECTS  ON  THE  PERIOD  IN  THE  ABSENCE  OF  GRAVITY 
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FIGURE  3-8.  COMPUTATED  EFFECTS  OF  THE  FREE  SURFACE  ON  THE  MAXIMUM 
EQUIVALENT  BUBBLE  RADIUS  IN  THE  ABSENCE  OF  GRAVITY 
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(B)  FINE  GRID  COMPUTATIONS 
FIGURE  3-9.  SPARK  BUBBLE  AT  SCALED  DEPTH  C  =  0.56 
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(A)  MEDIUM  GRID 
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FIGURE  3-11.  COMPUTED  AND  MEASURED  PLUME  HEIGHTS  FOR  THE 
C  =  0.56  SPARK  BUBBLE 
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CHAPTER  4 

COMPARISONS  TO  EXPERIMENTS 


In  this  chapter,  data  from  shallow  depth  explosion  experiments  are  presented  and  compared. 
These  tests  were  conducted  to  both  provide  additional  benchmark  validations  for  our  computational 
code  as  well  as  to  demonstrate  feasibility  for  creating  an  effective  water  barrier.37  Three  sets  of 
experiments  were  conducted  at  different  scales.  The  first  set  of  experiments  consisted  of  small 
(between  1-  and  2-lb)  individual  charges  detonated  at  the  NSWCCD  Test  Pond  in  May  1993.  Both 
above  surface  and  below  surface  photographs  were  taken  of  these  tests.  These  tests  will  be  denoted 
with  the  prefix  CD.  The  second  set  of  tests  were  performed  in  the  Arvonia  test  facility  (ARV)  and 
consisted  of  both  single  and  multiple  10-lb  charges  of  Composition  C-4.  The  experimental  results 
were  recorded  using  above  surface  cameras,  in  addition  to  the  use  of  conductivity  probes38  and 
Microwave  measurements.39  Larger  (25  lb)  charges  of  PBXN-103  were  tested  at  the  Aberdeen  Prov¬ 
ing  Grounds  Briar  Point  facility  (BP)  in  June,  1994.40,41  The  particular  shots  focused  on  are  listed  in 
Tables  4-1  and  4-2  below. 


TABLE  4-1.  SINGLE  CHARGE  TEST  SHOTS 


Shot 

Label 

Charge 

Type 

d 

(ft) 

W 

m 

4  max(°°) 

(ft) 

C 

CD-9 

C-4 

2.5 

1.25 

4.87 

0.51 

CD- 10 

Pentolite 

2.5 

1.72 

5.38 

0.46 

CD- 13 

Pentolite 

5.0 

1.73 

5.27 

0.95 

ARV-1 

C-4 

10.0 

6.0 

9.44 

0.64 

ARV-2 

C-4 

10.0 

10.0 

9.14 

1.09 

BP-8 

PBXN-103 

10.08 

26.3 

14.85 

0.67 

BP-9 

PBXN-103 

13.83 

26.3 

14.43 

0.96 

TABLE  4-2.  MULTIPLE  CHARGE  TEST  SHOTS 


Shot 

Label 

Charge 

Type 

d 

(ft) 

W 

(lb) 

Nc 

S 

(ft) 

M 

(Iblft) 

4  <2°) 
"me 

(ft) 

C 

ARV-3 

C-4 

8.0 

10.0 

5 

8.0 

1.25 

11.55 

0.69 

BP-12 

PBXN-103 

10.1 

26.3 

5 

10.1 

2.61 

20.83 

0.48 

BP-14 

PBXN-103 

13.8 

26.3 

5 

13.8 

1.90 

17.02 

0.81 

The  values  listed  for  the  scaled  depth  C  in  Table  4-1  were  computed  using  the  bubble  parameters 
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listed  in  Table  3-7,  together  with  equations  (3-7)  and  (3-14).  In  Table  4-2,  C  =  dIA where 
was  computed  using  (3 -18a)  using  a  cylindrical  bubble  assumption.  Also  listed  in  Table  4-2 
are  Nc ,  the  number  of  charges  used  for  the  shot,  S,  the  standoff  (distance)  between  charges,  and 
M  =  W/S,  the  mass  per  unit  length. 

COMPUTATIONAL  GRIDS 

The  grids  used  in  this  study  have  essentially  the  same  structure  as  those  used  in  the  spark  bub¬ 
ble  comparisons  described  in  Chapter  3.  That  is,  they  are  comprised  of  a  uniform  grid  region  of 
square  cells  that  extend  slightly  past  the  bubble  at  its  largest  size.  Outside  this  region,  the  grids  are 
extended  by  stretching  so  that  the  computational  boundaries  correspond  to  the  experimental  condi¬ 
tions,  whenever  possible.  For  example,  in  all  cases,  the  bottom  computational  boundary  corresponds 
to  the  location  of  the  bottom  below  the  charge.  Although  the  other  boundaries  of  the  test  sites  were 
not  axisymmetric  or  two-dimensional  they  were  generally  sufficiently  far  away  from  the  charge  as  to 
have  a  negligible  effect  on  the  bubble  dynamics. 

It  follows  from  (3-18)  that  the  maximum  to  minimum  bubble  radii  is  much  larger  in  the  two 
dimensional  cylindrical  case  than  for  spherical  bubbles.  For  example,  for  the  cases  listed  in 
Table  4-2,  the  values  of  a/2D)  corresponding  to  two-dimensional  approximations  to  shots  ARV-3, 
BP-12,  and  BP-14,  are  58.5,  44.1,  and  42.3,  respectively.  Using  a  single  uniform  grid  capable  of 
resolving  both  the  bubble  at  its  minimum  size,  while  extending  past  the  bubble  at  its  maximum  size 
would  require  an  excessive  number  of  grid  points.  However,  this  problem  can  be  alleviated  using 
two  separate  grids. 

Initially  a  grid  that  is  fine  in  a  region  surrounding  the  charge  was  used  until  the  bubble 
approached  the  boundary  of  the  fine  region.  Then  the  solution  was  remapped,  conserving  mass  and 
momentum,  onto  a  grid  that  was  coarser  than  the  fine  region  of  the  first  grid,  but  still  able  to  resolve 
the  bubble  after  the  initial  grid  had  been  used.  This  second  grid  was  uniform  in  a  large  enough 
region  to  contain  the  important  long  time  dynamics  of  the  problem.  For  example,  the  initial  grid  for 
shot  BP- 12  consisted  of  40  x  98  square  cells  with  size  h  =0.1  in  the  region  0  <  r  <  4.0  intersected 
with  -15.0  <  z  <  -5.2.  Outside  this  uniform  region  the  cells  were  stretched  radially  to 
r  =  Rl  =  160  using  40  additional  grid  lines,  and  downward  to  z  =  Zb  =  -54  using  30  grid  lines. 
Above  the  uniform  region  the  spacing  in  the  z  direction  was  uniformly  set  to  A z  =  0.2  in  the  region 

-5.2  <z<  4.0.  The  initial  grid  was  used  for  0  <  r  <  0.015  sec,  while  the  cylindrical  bubble  grew 

from  its  initial  radius  of  0.47  ft,  to  3.3  ft.  The  computed  solution  at  t  =  0.015  was  then  remapped 
onto  a  second  grid  having  square  cells  of  size  h  =  0.3  in  the  region  0  <  r  <24,  intersected  with 
-30  <  z  <30.  As  with  the  initial  grid  domain,  the  second  grid  was  stretched  radially  to 
r  =  RL  =160  using  40  grid  lines,  and  downwards  to  z  =  Zb  =  -54  using  30  grid  lines.  Above 
z  =  30  the  grid  was  also  stretched  upwards  using  the  stretching  A Zj  =  l.lAz/_ls  with  the  additional 
restriction  that  A  Zj  <15.  Using  50  grid  lines  above  the  uniform  region  allowed  the  computational 
domain  to  extend  to  z  =  ZT  =  326. 

While  the  initial  grid  is  used  for  only  a  small  fraction  of  the  first  bubble  period,  the  computa¬ 
tions  were  carried  out  for  several  periods  on  the  second  grid.  Typically,  the  use  of  the  initial  grid 
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represents  approximately  2  percent  of  the  total  computational  effort.  For  this  reason,  the  second 
grid  is  referred  to  as  die  “principle”  grid.  Since  it  was  also  observed  to  slightly  improve  die  single 
shot  computations  (it  typically  reduced  the  error  by  0.5  percent)  we  adopted  this  methodology  for 
all  the  runs.  A  summary  of  the  cell  sizes  in  the  uniform  regions  and  total  numbers  of  cells  (NXM) 
used  for  all  the  computations  are  listed  below  in  Table  4-3. 

TABLE  4-3.  SUMMARY  OF  COMPUTATIONAL  GRIDS  AND  RESULTS 


Bubble 

Maximum 

Radius 

Shot 

Initial  Grid 

Principle  Grid 

Period  (sec) 

Radius  (ft) 

Ratio 

Label 

h 

NxM 

h 

NxM 

T 

x  meas 

Th 

a  meas 
AMax 

Ah 

■^max 

a 

CD-9 

0.05 

70  X  150 

0.1 

80  X  275 

- 

0.180 

4.57 

4.53 

15.8 

CD- 10 

0.05 

70  X  150 

0.1 

80  X  280 

- 

- 

5.50 

5.44* 

16.5 

CD- 13 

0.05 

70  X  180 

0.1 

82  X  281 

0.212 

0.216 

5.75 

5.12 

16.1 

ARV-1 

0.1 

70  X  160 

0.2 

80  X  270 

0.34 

0.344 

- 

9.03 

15.3 

ARV-2 

0.1 

90  X  270 

0.2 

80  X  290 

- 

0.363 

- 

8.93 

14.8 

BP-8 

0.125 

74  X  175 

0.25 

96  X  299 

0.542 

0.540 

- 

14.47 

12.5 

BP-9 

0.125 

70  X  220 

0.25 

96  X  284 

0.544 

0.558 

- 

14.31 

12.1 

ARV-3 

0.05 

100  X  280 

0.2 

100  X  290 

- 

0.558 

- 

10.94 

58.5 

BP-12 

0.1 

80  X  174 

0.3 

120  X  280 

- 

0.872 

- 

19.30 

44.1 

BP-14 

0.1 

80  X  183 

0.3 

120  X  280 

- 

0.897 

- 

16.58 

42.3 

The  cell  sizes  in  the  principle  grid  uniform  region  were  selected  so  that  Amaxlh  ~  50, 
corresponding  to  the  “medium”  grid  used  for  the  spark  bubble  comparison,  and  the  h  =  0.2  grid 
listed  in  Table  3-1.  Therefore,  when  the  ratio  a  =  Amax/Amin  =  10,  it  can  be  expected  that  the  error 
in  computing  the  period  or  maximum  radius  will  be  approximately  2  percent.  For  the  single  shot 
cases  listed  in  Table  4-3,  the  values  of  a  are  slightly  larger  then  10,  and  the  error  can  be  expected 
to  be  approximately  2.5  percent.  The  initial  grid  uniform  cell  sizes  were  generally  selected  so  that 
Am\Jh  >  4  to  adequately  resolve  the  bubble  at  its  smallest  size. 

Also  listed  in  Table  4-3  are  the  computed  and  available  measurements  of  bubble  periods  and 
maximum  radii.  The  agreement  in  all  cases,  except  for  the  maximum  radius  for  shot  CD-13,  pro¬ 
vides  an  additional  justification  for  using  the  empirical  values  listed  in  Table  3-7.  With  the  excep¬ 
tion  of  shot  CD- 10,  which  is  expected  to  vent,  all  computations  were  performed  using  =  0.5  for 
0  <  t  <  0.8 Tdtb,  where  Tdb  was  computed  using  (3-19)  as  an  approximation  to  the  period.  For 
t  >  0.8 Td  b,  the  value  was  decreased  to  ea  =0.1  (see  Chapter  2).  For  shot  CD-10,  the  value  of  eA 
was  varied,  and  the  results  presented  in  the  next  section. 


*  Maximum  width  of  vented  cavity. 
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BUBBLE  PROFILE  COMPARISONS 

We  begin  our  validations  by  comparing  computational  bubble  shapes  and  sizes  to  the  underwa¬ 
ter  photography  performed  at  Carderock.  For  all  figures  in  this  chapter  the  units  of  distance  are 
measured  in  feet,  and  time  is  measured  in  seconds.  Figure  4-1  shows  the  computed  and  measured 
outlines  of  the  bubble  and  plume  at  t  =  0.1  sec  for  shot  CD-9,  which  roughly  corresponds  to  the 
time  of  maximum  bubble  size.  For  the  computations,  two  contours  of  density  are  shown,  namely 
p  =  0.9p0  and  p  =  0.0  lp0.  Note  that  the  lower  valued  contour  provides  a  better  match  for  the 
plume  outline.  The  fact  that  these  contours  are  widely  separated  in  the  plume  indicate  an  instability 
at  the  air-plume  interface.  This  is  because  of  the  plume’s  downward  acceleration  (the  bubble  has  a 
lower  pressure  than  the  air)  and  its  upward  velocity.  Inside  the  bubble;  however,  the  contours  agree 
closely  with  each  other  and  with  the  measurements.  The  quality  of  the  photographs  were  not  good 
enough  to  observe  the  downward  moving  jet  which  is  predicted  by  the  computations.  For  this 
example  A  =  4.87,  and  the  scaled  depth  C  =  0.51  is  close  to  the  critical  scaled  depth  when  vent¬ 
ing  can  be  expected.  The  computed  equivalent  maximum  radius  (based  on  maximum  volume)  was 
A  max  (C )  =  4.53.  According  to  (3-36)  the  expected  maximum  radius,  taking  the  free  surface  into 
account,  is  Amax(C)  =  4.61.  Using  the  same  grid  resolution  the  computed  maximum  radius  of  a 
spherically  symmetric  bubble  with  the  same  initial  conditions  (and  pressure  at  infinity  set  to  the 
hydrostatic  pressure  at  the  charge  depth)  was  A^ax  =  4.75,  representing  an  error  of  approximately 
2.5  percent.  This  explains  why  the  computed  bubble  is  slightly  smaller  than  the  measurements  and 
is  slightly  below  the  expected  value  of  Amax(C). 

An  even  shallower  case  was  shot  CD- 10  (W  =  1.72  lbs  of  Pentolite  at  a  depth  of  d  -  2.5  ft) 
which  has  a  scaled  depth  C  =  0.46.  In  this  case  the  bubble  is  expected  to  vent.  Figure  4-2  com¬ 
pares  the  computed  and  measured  bubble  outlines  at  t  =  0.10  sec  using  5  computations  with 
different  values  of  the  air-water  interface  cutoff  eA  (cf.  (2-3b)).  When  zA  =  0.5  (Figure  4-2(A)) 
venting  did  not  occur  and  the  computed  bubble  is  substantially  smaller  than  the  observations.  Using 
eA  -  0-4  causes  the  bubble  to  eventually  vent  at  time  Tv  =  0.076,  but  the  bubble  profile  is  not 
greatly  effected  at  t  =0.10  (Figure  4-2(B)).  However,  the  central  water  plume  is  slightly  narrower 
and  the  density  contour  p  =  0.01p0  is  higher  than  in  the  previous  case.  As  eA  is  decreased  (Figure 
4-2(C-E))  the  bubble  vents  earlier,  causing  the  cavity  to  get  larger.  In  Figures  4-2  (D)  and  (E)  the 
outline  of  the  cavity  is  in  good  agreement  to  the  photographic  observations.  Also,  the  amount  of 
water  in  the  central  region  is  greatly  decreased  as  the  bubble  vents  earlier.  However,  as  mentioned 
in  Chapter  2  the  venting  of  the  bubble  is  crudely  approximated  in  our  model,  with  the  pressure  in 
the  cavity  instantaneously  changing  from  the  value  it  had  before  the  bubble  vents  to  the  constant 
and  uniform  atmospheric  pressure.  Therefore,  long  time  predictions  of  bubbles  which  have  vented 
cannot  be  expected  to  be  reliable. 

PLUME  COMPARISONS 

At  early  times  after  detonation  the  plume  outlines  were  observed  clearly  in  the  photographic 
records.  Comparisons  of  the  observed  plume  heights  and  widths  provide  a  critical  validation  of  the 
computational  model. 
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Single  Charge  Plume  Heights 

Figure  4-3  shows  the  measured  and  computed  values  of  the  plume  heights  for  the  Carderock 
shots  CD-9,  CD- 10,  and  CD-13.  The  computations  for  CD-9  were  the  same  as  shown  in  Figure  4-1. 
The  plume  heights  for  CD- 10  were  taken  from  the  same  computation  depicted  in  Figure  4-2(C) 
(with  eA  =  0.3).  For  the  computations,  the  plume  height  was  defined  to  be  the  highest  z  value  for 
which  p  >  0.01p0.  The  “ripples”  in  file  computed  plume  heights  are  due  to  the  coarse  stretched 
grids  which  were  used  above  z  =  13  ft.  While  the  agreement  for  the  shallow  cases  (CD-9  and 
CD-10)  is  very  good,  the  computation  for  shot  CD-13,  where  the  scaled  depth  is  C  =  0.95,  underes¬ 
timates  the  measured  plume  heights  by  over  30  percent. 

The  underestimation  of  plume  heights  (using  an  incompressible  liquid  assumption)  has  been 
discussed  by  Kedrinski,5  who  also  noted  that  early  time  plumes  will  be  enhanced  if  the  surface  is 
indented  above  the  charge.  Such  a  recess  would  be  caused  by  the  spallation  region  caused  by  the 
interaction  of  the  shock  with  the  free  surface.  Since  our  model  uses  an  incompressible  assumption 
for  water,  this  phenomenon  cannot  be  predicted  directly  using  our  methodology.  However,  we  can 
simulate  this  shock  effect  by  modifying  the  initial  conditions.  For  example,  instead  of  initializing 
the  free  surface  to  be  flat  (z  =  0),  we  can  recess  the  surface  by  “cutting  out”  a  spherical  region 
with  radius  i?;  centered  at  z  =d,  where  z  =-d  is  the  initial  charge  location.  The  computed  plume 
heights  for  shot  CD-13,  repeated  with  an  initial  recess  radius  =  1.03 d,  are  also  displayed  in  Fig¬ 
ure  4-3.  The  discrepancy  between  the  measured  and  computed  data  in  this  case  is  greatly  reduced. 

Computational  grids  used  for  the  Arvonia  shots  were  twice  the  size  as  those  used  for  the  Car¬ 
derock  computations  so  that  the  relative  resolution  A^/h  was  the  same  as  for  shot  CD-9. 
Figure  4-4  shows  comparisons  of  measured  and  computed  plume  heights  for  shots  ARV-1  and 
ARV-2.  In  each  case  the  computations  with  an  initially  flat  free  surface  underpredicted  the  plume 
heights  by  30  percent  to  as  much  as  50  percent.  As  with  shot  CD-13,  this  discrepancy  was  greatly 
reduced  with  the  choice  Rt  =  1.03 d,  for  both  cases. 

Figure  4-5  compares  the  initial  plume  formation  for  shot  ARV-1  based  on  computations  using 
a  recessed  and  flat  initial  free  surfaces.  Here,  the  density  contour  p  =  p0  is  displayed.  In  the 
recessed  case,  a  thin  jet  of  water  forms  at  an  early  time  and  precedes  the  thicker  “main”  column 
which  would  appear  without  the  recess.  This  recess  has  little  or  no  influence  on  the  bubble  dynam¬ 
ics  under  the  surface. 

The  empirical  value  R(  =  1.03<i  was  also  used  for  charges  of  PBXN-103.  Figure  4-6  includes 
measured  and  computed  plume  heights  for  shots  BP-8  and  BP-9  performed  at  the  Briar  Point  Facil¬ 
ity  of  the  Aberdeen  Proving  Grounds.  These  shots  were  conducted  with  25  lb  charges  of 
PBXN-103  with  a  2.5  lb  C-4  booster.  This  is  roughly  equivalent  to  a  26.3  lb  charge  of  PBXN-103. 
For  shot  BP-9,  where  C  =  0.96,  the  plume  heights  are  substantially  underestimated  with  Rt  =  0.  For 
the  early  times,  0  <  t  <  1.0,  the  computed  heights  with  Rt  =  1.03  are  20  to  25  percent  below  the 
measured  values.  However,  this  discrepancy  disappears  at  later  times.  For  shot  BP-8  (C  =  0.67)  the 
initially  flat  surface  yielded  predictions  that  were  over  50  percent  under  the  measured  values.  As 
with  shot  BP-9,  the  computations  using  Rt  =  1.03 d  underestimated  the  measured  values  at  the  early 
times  (0  <  t  <  0.7).  We  conjecture  that  this  underestimation  is  caused  by  the  initial  shock  dome  of 


4-5 


NSWCDD/TR-94/156 


spalled  water  which  is  not  modeled  in  our  simulations.  At  later  times  the  computed  heights  become 
slightly  greater  than  the  measured  values.  Other  effects,  such  as  the  plume  breakup  and  drag  forces, 
may  become  significant  at  these  times. 

Based  on  these  results  as  well  as  those  displayed  in  Figures  3-11  and  4-3,  the  value 
Rt  ~  1.03 d  appears  to  provide  a  reasonably  good  match  to  plume  height  data  for  scaled  depths  C  in 
the  approximate  range  0.6  <  C  <  1.1  for  both  C-4,  Pentolite,  and  PBXN-103  charges.  The  use  of 
the  initially  indented  surfaces  had  a  negligible  effect  on  the  computed  periods  and  maximum  bubble 
radii  (<  0.1  percent)  for  the  cases  considered.  The  computed  values  listed  in  Table  4-3  were 
unaffected  by  the  use  of  the  indented  initial  surfaces. 


Multiple  Charge  Plume  Heights 

We  now  turn  our  attention  to  the  multiple  charge  shots  listed  in  Table  4-2.  In  the  cases  con¬ 
sidered  the  distance  between  the  individual  charges  was  less  than  the  free  field  maximum  radius  of 
the  single  charges.  For  example,  a  single  10-lb  charge  of  C-4  at  a  depth  of  8  ft  has  a  free  field 
maximum  bubble  radius  of  9.29  ft,  which  is  larger  than  the  standoff  distance  between  charges  of 
8  ft  for  shot  ARV-3.  The  corresponding  single  charge  values  for  shot  BP-12  and  BP-14  are  listed 
in  Table  4-1  (Amax(oo)  values  for  shots  BP-8,  and  BP-9,  respectively).  In  such  cases,  it  can  be 
expected  that  the  individual  bubbles  will  merge  into  a  nearly  cylindrical  bubble  (given  a  sufficient 
number  of  charges).  Therefore,  a  two-dimensional  approximation  may  be  reasonable  near  the  center 
of  the  line  of  charges. 

Figure  4-7  compares  the  measured  and  computed  plume  heights  for  shots  BP-12,  and  BP-14. 
For  a  simultaneous  detonation  of  multiple  charges,  spallation  can  be  expected  not  only  above  each 
charge,  as  in  the  single  shot  case,  but  also  in  between  the  charges  due  to  the  interaction  of  the  shock 
intersections  with  the  free  surface.  Shock  interactions  from  underwater  spherical  charges  have  been 
studied  by  Stebnovskii42  and  by  Colebum  and  Roslund,43  but  not  in  conjunction  with  the  air-water 
surface.  This  interaction  effect  is  clearly  seen  in  the  images  taken  from  video  tapes  of  shots  BP-12 
and  BP-14.  A  sequence  of  these  frames  appear  in  the  top  half  of  Figure  4-8  (shot  BP-12)  and  Fig¬ 
ure  4-9  (shot  BP-14).  In  each  case  four  distinct  “peaks”  are  clearly  observed  at  early  times,  which 
jet  upwards  from  the  midpoints  of  the  five  underwater  charges.  The  fact  that  this  effect  can  be  seen 
in  the  first  video  frame  after  the  detonation  (at  t  =  0.03)  supports  the  conjecture  that  this  is  a  shock 
effect.  In  Figure  4-7,  both  the  measured  heights  of  the  peaks,  and  the  “valleys”  are  plotted  for 
shot  BP-14.  The  computations  are  closer  to,  but  still  below  the  measured  heights  of  the  valleys. 
This  discrepancy  is  probably  caused  by  the  spallation  directly  above  each  charge  as  occurred  with 
the  single  shot  cases  previously  discussed. 

As  with  the  single  shot  cases,  multiple  shot  shock  effects  can  be  modeled  empirically  by  ini¬ 
tializing  the  air  water  surface  with  recesses  corresponding  to  where  spallation  occurs.  In  addition  to 
a  recess  directly  above  each  charge,  the  surface  should  also  include  recesses  in  between  the  charges, 
where  shock  focusing  is  expected.  This  latter  recess  could  be  modeled  as  the  region  of  intersection 
between  spheres  centered  at  reflected  points  located  a  distance  2d  above  each  charge.  This  model 
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is  shown  in  Figure  4-10.  Testing  this  model  computationally  requires  the  use  of  our  three- 
dimensional  code  BUB3D  which  will  be  done  in  a  forthcoming  paper. 


Secondary  Plumes  -  Long  Time  Behavior 

To  better  understand  the  quantitative  comparisons  of  plume  densities,  we  must  also  qualita¬ 
tively  examine  the  long  time  plume  dynamics.  The  long  time  plume  behavior  can  be  seen  in  Figures 
4-8  and  4-9  for  multiple  shots  BP-12,  and  BP-14,  and  in  Figures  4-11,  and  4-12  for  single  shots 
BP-8  and  BP-9.  The  computational  results  for  all  four  cases  are  shown  using  gray  scale  density 
contours.  A  middle  value  of  gray  was  used  to  shade  the  region  0.5p0  <  p  <  p0.  For  values  of 
p  <  0.5p0  the  gray  scale  was  lightened  using  linear  interpolation,  so  that  p  =  0  corresponds  to  white. 
The  size  scale  used  for  the  contour  frames  is  roughly  the  same  as  in  the  photographic  frames  in 
each  case.  As  indicated  in  the  computational  frames,  the  top  of  the  video  frames  corresponds  to  a 
height  of  approximately  90  ft.  We  remark  that  in  regard  to  the  multiple  shot  cases  in  Figures  4-8 
and  4-9,  the  photographic  view  is  at  nearly  a  right  angle  to  the  computational  contour  of  the  plume 
and  bubble  cross  section. 

The  single  shot  examples  shown  in  Figures  4-11  and  4-12  begin  similarly  as  the  previous 
examples  shown  in  Figures  3-9,  4-1,  and  4-5.  Here,  two  central  jets,  a  large  one  rising  into  the 
atmosphere  and  a  smaller  one  falling  into  the  bubble,  form  between  the  bubble  and  air  surface  as 
the  bubble  begins  its  first  contraction.  The  downward  moving  jot  can  be  observed  in  the  computed 
frame  labeled  T  -  0.41  in  Figure  4-11.  The  downward  jet  is  just  beginning  to  form  at  time 
T  =  0.41  for  shot  BP-9  (Figure  4-12)  and  by  T  =  0.61  has  collided  with  another  thicker  jet  moving 
upward  from  the  bottom  of  the  bubble.  For  a  further  discussion  on  the  formation  of  jets  in  bubbles 
near  a  free  surface  see  Blake  et.  al.,34  who  also  provides  an  analysis  of  such  phenomenon  through 
the  use  of  the  Kelvin  impulse  concept. 

The  photographs  of  single  shot  BP-8  in  Figure  4-11  show  a  persistent  central  column,  which 
appears  to  be  initially  surrounded  by  spray,  but  at  later  times  appears  to  be  a  thin  column  of  water. 
Although  this  column  appears  to  be  thinner  in  the  computations,  its  height  and  persistence  are 
matched  very  well,  even  through  the  last  frame  at  time  t  =  5.8  sec.  The  long  time  plume  heights 
are  also  reproduced  accurately  by  the  computations  for  shot  BP-9,  as  shown  in  Figure  4-12. 

An  interesting  effect  that  appears  both  in  the  computation  and  photographs  is  the  emergence  of 
secondary  plumes  ejected  upward  and  radially  away  from  the  central  column.  These  plumes  can  be 
seen  emerging  during  the  bubble’s  second  expansion.  In  Figure  4-11,  at  time  t  ~  0.8  the  ejection  of 
these  plumes  can  be  seen  in  both  the  computation  and  the  photograph.  The  computed  and  photo¬ 
graphed  heights  of  these  secondary  plumes  are  in  good  agreement,  although  the  widths,  (or  angles 
from  the  vertical)  are  somewhat  narrower  in  the  computation.  The  secondary  plumes  appear  a  little 
later  in  the  photographs  of  shot  BP-9  (Figure  4-12).  In  this  case,  the  fallback  of  the  radial  plumes 
were  matched  well  by  the  computation,  but  the  width  and  height  are  slightly  overestimated.  The 
ability  to  capture  the  secondary  plume  phenomena,  even  qualitatively,  is  a  difficult  task  due  to  the 
complex  structure  of  the  bubble  surface  and  the  relatively  long  time  it  takes  to  develop. 
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Note  that  these  secondary  plumes  also  appear  in  the  multiple  charge  shots  displayed  in  Figures 
4-8  and  4-9.  The  emergence  of  these  plumes  is  evident  by  time  t  =  1.2  for  both  shots.  The  two 
dimensional  computations  predict  the  emergence  of  the  secondary  plumes  at  times  roughly  coincid¬ 
ing  with  the  photographs.  For  both  of  these  multiple  charge  shots,  the  central  plume  falls  faster  in 
the  computation  than  in  the  experiment,  particularly  for  shot  BP- 14.  This  can  be  expected  since  the 
central  plume  heights  are  underpredicted.  As  mentioned  earlier,  we  expect  this  discrepancy  to  be 
reduced  by  including  the  three-dimensional  shock  effects,  using  the  empirical  model  described  in 
Figure  4-10.  Another  phenomenon  that  will  effect  the  duration  of  the  plume  is  the  drag  on  the 
small  water  droplets  as  they  fall  with  gravity.  As  the  plume  breaks  up  into  small  droplets  or  spray, 
this  effect  can  become  significant.  Droplet  breakup  and  drag  effects  are  currently  not  included  in 
our  model. 


MICROWAVE  COMPARISONS 

Quantitative  measurements  of  the  plume  density  were  made  using  microwave  measurements. 
These  measurements  were  based  on  the  amount  of  microwave  absorption  through  the  plume. 
Microwaves  were  sent  and  received  using  a  pair  of  three  foot  radius  parabolic  dishes  placed  on 
either  side  of  the  plume  at  equal  heights  above  the  water  surface.38  To  compare  the  microwave 
measurement  with  the  computed  densities,  the  computed  values  were  integrated  within  a  cylindrical 
region  of  radius  3  ft  between  the  microwave  dishes.  In  Cartesian  (x,y,z)  coordinates  this  integral 
is 

Hq+Rd  R(z)  oo 

In  =  J  J  Jp n(x,y,z)  dx  dy  dz  (4-1) 

Hd  —Rj)  —R  (z  )  — 

where  HD  is  the  height  of  the  center  of  the  dish,  RD  is  the  radius  of  the  dish,  pn  is  the  density  at 
time  tn ,  and 

U&D  -  (z-Hd)2  if  \z  —  Hd  I  <  Rd , 

~\  0  otherwise.  ($-2) 

The  integral  In  corresponds  to  the  total  mass  of  water  in  the  cylindrical  region  between  the  two 
dishes  straddling  the  plume  at  time  tn . 

For  the  two-dimensional  approximations,  p”  does  not  change  in  the  y  direction  (parallel  to  the 
line  of  charges,  and  perpendicular  to  the  line  between  the  microwave  dishes).  By  symmetry  across 
x  =  0  it  follows  that 

Hd+Rd  CO 

In  =  4  J  R(z)jpn(x,z)  dx  dz.  (4-3) 

Hd-Rd  0 

This  integral  is  approximated  using  the  midpoint  rule  quadrature 

Iu  =  (Ac),  (A Z)j,  (4-4) 

j=J  1  i'=l 
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where  p  ”j  is  the  computed  cell  centered  density  of  cell  (*,  j), 


zj  = 


zf  +  zf 


2  ,  (&z)j  =  zj  -  zf,  (Axf  =xi+ j  -  xt , 

zj  -  ram(HD  +  RD ,  zj+i\  zf  =  max(i7D  -  RD ,  zf), 

and  Ji  and  J2  are  indices  such  that  zJl  <  HD  -  RD ,  and  Zj2  >HD  +RD.  In  the  above,  cell  (i,j) 
corresponds  to  the  rectangular  region  xt  <x  <  xi+1,  Zj  <  z  <  zj+l  for  i  =  1, N  and  j  =  1, M. 
In  cylindrical  (r,  0,  z)  coordinates  the  integral  (4-1)  can  be  written  as 


Hd+Rd  2  R(z)lcos(Q) 

/”  =  4  J  J  J  r  pn(r,z)dr  dQ  dz. 


(4-5) 


Hd-Rd  0=0  r=0 


This  integral  is  approximated  on  the  computational  grid  using  the  quadrature 


J 2  Afe 

4"  =  4£  £ 

j=J i  i=i 


(A0),  mj. 


(4-6a) 


where 


ri  = 


rP  +  rf 


2 

,R- 


(Ar),.  =  rf  -  rt,  &,=(/-  (49),  = 


7t 


2 Me 


2Mc 


rf  =  min(R;  / ,  ri+1),  rf  =  min(i?jV ,  r{ ),  Rjj  = 


R(Zj) 


(4-6b) 


(4-6c) 


COS(0j) 

Note,  that  the  values  ft  depend  implicitly  on  the  indices  l  and  k  as  indicated  by  (4-6b)  and  (4-6c) 
due  to  the  variable  limits  of  integration  in  the  r  coordinate. 

For  the  actual  comparisons,  the  total  mass  values  In  were  scaled  by  both  the  cross  sectional 
area  of  the  cylinder  and  the  water  density  to  obtain  an  “effective  water  length”  (EWL).  This  value 

In 


EWLn  = 


Po  t&d 


(4-7) 


corresponds  to  the  length  of  water  filling  the  cylinder  having  an  equivalent  mass  as  the  plume  at 
time  tn . 


Before  proceeding  with  the  comparisons,  note  that  because  of  the  novelty  of  using  microwave 
absorption  to  measure  plume  densities,  no  independent  validation  of  this  procedure  is  currently 
available.  Therefore,  these  comparisons  should  be  considered  in  conjunction  with  the  plume  height 
data,  photographs,  measured  bubble  periods,  and  the  probe  measurements  in  the  following  section, 
whenever  these  were  available.  In  this  way,  the  microwave  comparisons  provide  validations  not  only 
for  the  computations,  but  for  the  measurements  as  well. 
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Figures  4-13  and  4-14  display  measured  and  computed  microwave  data  for  the  single  charge 
shots  ARV-1  and  ARV-2.  For  all  the  Arvonia  shots,  the  microwave  measurements  were  taken  at  a 
height  Hd  =  25  ft.  The  measured  data  for  shot  ARV-1  show  an  initial  rise  in  the  EWL  from  0  to 
0.2  occurring  between  0.2  and  0.25  sec  after  detonation.  There  is  an  inflection  point  in  the  meas¬ 
ured  data  at  approximately  t  =  0.25,  after  which  the  EWL  spikes  to  a  value  of  approximately  1.6  at 
time  t  =  0.55.  Using  an  initially  flat  surface,  the  computed  data  (Rt  =  0)  misses  the  inflection  point 
entirely,  and  the  EWL  increases  from  0  to  approximately  0.8  in  the  time  interval  0.3  <  t  <  0.6. 
However,  the  computed  EWL  using  an  initially  recessed  surface  (R(  =  1.03c?)  indicates  an  initial 
rise  from  0  to  0.1  during  the  time  interval  0.2  <  t  <  0.3  with  an  inflection  point  at  approximately 
t  =  0.35.  Based  on  this  data,  it  is  reasonable  to  expect  that  the  first  inflection  point  (or  kink)  in  the 
measured  data  is  due  to  shock  effects  (cf.  Figures  4-4  and  4-5).  Furthermore,  the  empirical  model 
used  to  simulate  these  effects  appears  to  underestimate  the  EWL  before  the  initial  kink.  This  can  be 
expected  since  the  spray  dome  that  is  spalled  upward  due  to  the  shock  reflection  is  not  included  in 
the  computations.  At  later  times  both  computations  show  a  peak  in  the  EWL  at  t  ~  1,  with  the 
R{  =  0  computation  predicting  EWL  values  approximately  30  percent  larger  than  the  R{  =  1.03 d 
computation  in  the  time  interval  1.0  <  t  <  2.2.  The  computed  peaks  in  the  EWL  occur  when  the 
secondary  plumes  pass  through  the  microwave  cylinder.  Even  though  this  event  was  not  repro¬ 
duced  in  the  measured  data,  the  overall  agreement  is  reasonable,  at  least  from  the  point  of  view  of 
giving  the  correct  order  of  magnitude  for  the  EWL  most  of  the  time.  At  late  times  the  measured 
data  lies  above  the  computed  data. 

Figure  4-14  displays  a  more  dramatic  difference  in  the  computations  using  a  flat  and  indented 
initial  free  surface.  Here,  the  flat  surface  computation  yields  completely  erroneous  results,  indicat¬ 
ing  no  plume  density  until  after  t  -  1.6  sec  while  the  measurements  and  observed  plume  heights  (cf. 
Figure  4-4)  indicate  that  the  plume  enters  the  microwave  cylinder  at  approximately  t  =  0.6.  Furth¬ 
ermore,  the  Rt  =  0  computation  also  predicts  a  rise  in  EWL  values  between  t  =  1.6  and  t  =  3.0. 
This  effect  was  caused  by  a  rebounding  jet  that  formed  after  the  computed  bubble  vented  shortly 
after  the  time  of  the  second  bubble  maximum.  However,  the  computed  EWL  values  with  the  ini¬ 
tially  recessed  surface  (Rt  =  1.03d)  are  in  excellent  agreement  with  the  measurements.  The  slight 
delay  (approximately  0.2  sec)  in  the  initial  rise  and  decline  of  computed  EWL  values  corresponds 
to  the  underestimation  of  initial  plume  heights  as  indicated  in  Figure  4-4.  The  difference  between 
the  measured  and  computed  peak  values  in  this  case  is  only  approximately  12  percent.  Also,  unlike 
the  Rt  =  0  computation,  the  bubble  did  not  vent  after  the  second  bubble  maximum,  and  no  rebound 
jet  formed.  This  example  demonstrates,  that  in  certain  cases,  particularly  for  values  of  C  =  1,  even 
the  long  time  dynamics  can  be  sensitive  to  relatively  minor  changes  in  the  initial  conditions  (e.g., 
recessed  surfaces  caused  by  shock  effects). 

Figure  4-15  compares  the  measured  and  computed  microwave  data  for  multiple  shot  ARV-3 
(cf.  Table  4-2).  The  microwave  equipment  was  calibrated  only  for  values  EWL  <  1.58  for  these 
experiments  as  indicated  by  the  plateau  in  the  measured  data  for  0.7  <  t  <  2.7.  Despite  the 
underprediction  by  the  computations  between  2.0  <  t  <  3.0,  once  again  the  agreement  is  good.  The 
initial  rise  times  agree  to  within  0.05  sec.  The  two  peaks  in  the  computations  at  approximately 
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t  =0.7,  and  t  =  1.35  correspond  to  the  rising  of  the  central  plume  followed  by  the  rise  of  the 
secondary  plumes,  through  the  microwave  cylinder.  The  “bumps”  at  the  later  times  correspond  to 
parts  of  the  plume  structure  falling  back  through  the  cylinder  toward  the  surface.  As  in  Figure  4-13, 
the  measured  EWL  values  take  longer  to  decay  to  0,  possibly  because  of  the  unmodeled  drag 
forces. 

Figure  4-16  compares  the  computed  and  measured  EWL  for  shot  BP-8  at  the  two  heights 
Hd  =  25  and  HD  =  12.5.  The  results  using  Rt  =  1.03 d  are  slightly  better,  primarily  at  early  times 
at  Hd  =  25.  At  this  height  at  later  times,  the  computations  indicate  that  there  is  approximately 
twice  as  much  water  as  measured  using  the  microwave  absorption  data.  This  is  somewhat  surprising, 
when  the  photographs  of  Figure  4-11  are  taken  into  consideration.  As  noted  earlier,  the  computed 
plumes  appear  to  be  substantially  thinner  than  in  the  photographs.  If  the  microwave  measurements 
are  accurate,  the  actual  plume  must  not  be  very  dense,  even  through  the  middle.  However,  the  per¬ 
sistence  of  the  central  plume  suggests  a  relatively  dense  core  in  the  plume,  as  the  computations 
predict.  It  is  also  possible  that  the  microwave  measurements  miss  part  of  the  plume  at  later  times 
due  to  wind  forces.  A  slight  drift  in  the  plume  is  noticeable  in  the  photographs,  particularly  during 
its  descent.  At  the  height  HD  =  12.5  ft,  the  microwave  measurements  saturated  at  EWL  =  0.9.  This 
comparison  only  indicates  reasonably  good  agreement  between  rise  and  fall  times  for  the  EWL 
values. 

Figure  4-17  compares  the  computed  and  measured  EWL  values  for  shot  BP-9.  For  this  case 
there  is  a  large  difference  between  the  R;  =  0  and  R{  =  1.03d  cases.  When  =  0  the  central  jet 
does  not  quite  reach  a  height  of  25  ft,  and  the  only  measured  water  is  produced  by  the  secondary 
plumes.  The  computation  with  2?f  =  1.03d  indicates  an  initial  rise  in  EWL  values  after  t  =  0.4, 
about  0.2  sec  later  than  the  microwave  data.  This  is  consistent  with  the  times  that  the  plume 
heights  approach  25  ft  as  indicated  in  Figure  4-6.  The  second  EWL  peak  at  t  =  1.0  when  HD  =  25 
is  due  to  the  radial  plume  passing  through  the  microwave  cylinder.  The  radial  plumes  shown  in  the 
photographs  of  Figure  4-12  do  not  appear  to  contribute  to  the  microwave  data,  and  perhaps  do  not 
reach  the  height  of  25  ft.  In  the  computations  the  radial  plumes  rise  above  30  ft  and  contribute 
significantly  to  the  EWL  values.  This  partly  explains  the  overestimation  of  the  computed  EWL 
values  for  t  >  1.2  sec.  However,  the  photographs  in  Figure  4-12  shows  a  substantial  amount  of 
water  in  the  plume  at  time  t  =  2.53  sec,  which  the  microwave  data  fails  to  indicate. 

At  Hd  =  12.5  both  of  the  computations  and  the  measurements  are  in  good  agreement  for  the 
initial  rise  of  EWL  values  between  0.2  <  t  <  0.4.  The  computations  of  the  EWL  values  diverge 
after  t  =  1.4.  This  difference  is  similar  to  what  occurred  with  the  analogous  computations  for 
ARV-2  where  the  bubble  vents  shortly  after  its  second  bubble  maximum  in  the  Rt  =0  case.  The 
high  EWL  values  in  this  case  are  caused  by  the  rebound  jet  passing  through  the  microwave  cylinder 
at  Hd  =  12.5.  When  Rt  =  1.03c?  the  venting  does  not  occur  and  the  agreement  with  the  measure¬ 
ments  is  much  better.  We  remark  that  there  is  no  evidence  from  the  photographs  of  this  rebound  jet 
(cf.  Figure  4-12).  However,  as  with  shot  BP-8,  the  measured  EWL  values  appear  to  fall  too 
quickly,  since  the  observed  plume  is  above  HD  =  12.5  until  some  time  after  t  =  3,  while  the  meas¬ 
ured  EWL  values  are  negligible  for  t  >  2.2. 
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Figure  4-18  compares  the  computed  and  measured  EWL  values  for  multiple  charge  shot  BP-12. 
At  both  heights  the  initial  rise  times  are  predicted  very  well,  while  the  computed  values  decay  ear¬ 
lier  (as  occurred  with  shot  ARV-3).  Since  the  early  time  EWL  values  agree,  while  the  computed 
plume  heights  were  below  the  measurements  (cf.  Figure  4-7),  it  can  be  expected  that  the  peaks  of 
the  plumes  at  early  times  do  not  contain  a  significant  amount  of  water.  Also  note  that  the 
microwave  dishes  were  centered  above  the  central  charge,  which  corresponds  to  the  central  valley  of 
the  initial  plume  structure  (cf.  Figure  4-8).  The  overall  agreement  between  the  computations  and 
measurements  in  this  case  is  good. 

Finally,  Figure  4-19  compares  the  computed  and  measured  EWL  values  for  multiple  shot 
BP-14.  Here  the  measured  EWL  rise  times  again  show  good  agreement  with  the  computed  values. 
The  computed  EWL  rise  appears  a  little  early  at  HD  =  12.5  and  a  little  late  at  HD  =  25.  The  local 
peaks  at  t  ~  1.2  are  caused  by  the  secondary  plumes  passing  through  the  microwave  cylinders  (cf. 
Figure  4-9).  The  peak  value  for  the  EWL  at  HD  =  12.5  of  22  ft  at  t  ~  0.5  is  clipped  off  in  the 
figure. 


DENSITY  PROBE  COMPARISONS 

A  second  set  of  plume  density  measurements  were  obtained  through  the  use  of  conductivity 
probes.  These  probes  were  first  developed  by  Phillips  and  Scott,43  and  consist  of  two  parallel  stain¬ 
less  steel  rods.  Phillips  and  Scott  found  that  the  conductivity  through  each  probes  was  linearly  pro¬ 
portional  to  the  unwetted  length  of  rod.  These  probes  were  originally  used  for  measuring  bubble 
radii  from  underwater  explosions. 

Lipton38  suspended  conductivity  probes  above  the  surface  to  obtain  plume  density  measure¬ 
ments  for  both  the  Arvonia  and  Briar  Point  tests.  Comparing  the  computed  density  at  a  specific 
location  to  the  measured  density  from  a  single  probe  has  little  meaning  for  at  least  two  reasons: 

(1)  The  free  surface,  and  hence  plumes,  undergo  periods  of  Rayleigh-Taylor  instability. 
Therefore,  point  measurements  cannot  be  expected  to  be  reproducible.  This  is  supported 
by  the  fact  that  probes  located  the  same  distance  but  on  opposite  sides  of  the  charge 
center,  yielded  vastly  different  density  histories. 

(2)  The  probes  were  suspended  by  a  rope  in  the  Arvonia  tests  which  was  deflected  several 
feet  upon  impacts  from  the  plumes.  Therefore,  the  actual  probe  locations  were  not 
known. 

Despite  the  instability,  some  degree  of  reproducibility  can  be  expected  from  an  integral  norm.  As 
with  the  microwave  comparisons,  the  EWL  value  can  be  used,  where  the  integration  over  the 
cylinder  (4-1)  is  replaced  by  integration  on  a  line.  That  is, 

co 

EWL "  =  J  p(*  ,yp  ,zp  )dx ,  (4-8) 

r0  — oo 

representing  integration  over  the  horizontal  line  located  at  y  =yp,  and  z  =  zp.  Note,  that  (4-8)  is 
simply  the  limit  of  (4-7)  and  (4-1)  as  RD  —>  0.  Since  the  probes  were  located  at  discrete  points,  a 
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trapezoid  rule  integration  was  used  to  determine  the  measured  EWL .  That  is 


EWL£  = 


n- 1 


Pj/2+  £p  i 
i=2 


(4-9) 


where  Sp  is  the  uniform  spacing  between  the  n  probes  and  p;  is  the  density  at  the  ith  probe  loca¬ 
tion. 


Comparisons  with  the  simulations  were  made  using  the  same  formula  (4-9),  except  the  density 
values  were  interpolated  from  the  computed  cell  density  values.  Using  (4-9)  on  the  computed  data 
at  the  same  locations  as  the  probes  is  referred  to  as  “Probe  Integration”  (P.I.).  In  addition,  the 
computed  values  were  integrated  over  the  same  grid  that  was  used  for  the  simulation.  This  “Full 
Integration”  (F.I.)  was  used  to  determine  if  there  was  plume  structure  missed  by  the  probes. 

For  the  Arvonia  tests  eight  probes  were  placed  on  a  rope  at  a  height  12.5  ft  above  the  surface. 
The  probes  were  placed  5  ft  apart.  Comparisons  between  the  measured  EWL  values  and  the  Com¬ 
puted  values  using  both  probe  and  full  integration,  are  displayed  in  Figures  4-20,  4-21,  and  4-22,  for 
shots  ARV-1,  ARV-2,  and  ARV-3,  respectively.* 

In  each  case,  the  computed  and  measured  initial  rise  times  for  the  EWL  values  are  in  good  agree¬ 
ment.  For  shots  ARV-1,  and  ARV-2,  this  can  be  expected  because  of  the  agreement  for  the  plume 
heights  displayed  in  Figure  4-4.  The  data  displayed  for  ARV-1  and  ARV-2  in  Figures  4-20  and  4- 
21,  correspond  to  the  computations  with  the  initially  recessed  surface  {Rt  =  1.03d).  Notice  the  F.I. 
values  rise  earlier  than  the  P.I.  values  for  these  two  cases.  This  is  due  to  the  fact  that  the  central 
plume  with  the  recessed  surface  is  thin  (see  Figure  4-5).  Therefore,  it  is  not  initially  detected  by  the 
density  values  at  the  probe  locations  which  are  2.5  ft  on  either  side  of  the  charge  center.  The  com¬ 
puted  and  measured  rise  time  of  EWL  values  for  multiple  shot  ARV-3  is  excellent.  This  suggests 
that  the  height  of  the  initial  plume  containing  substantial  water  mass  is  accurately  prediced  by  the 
computations. 

It  is  obvious  from  these  figures  that  there  is  little  agreement  for  the  EWL  values.  The  large 
discrepancy  between  the  computed  P.I.  and  F.I.  values  suggests  that  the  probe  spacing  is  much  too 
large  to  accurately  portray  the  plume  structure.  Furthermore,  it  is  also  unlikely  that  the  measured 
results  are  very  reproducible,  since  the  integration  used  for  the  EWL  usually  involved  contributions 
from  only  3  or  4  of  the  probes.  Finally,  the  high-frequency  oscillations  in  the  measured  data,  also 
make  comparisons  difficult,  unless  these  are  substantially  smoothed. 

For  the  Aberdeen  tests  the  number  of  probes  was  increased  to  16,  and  the  spacing  between  the 
probes  was  reduced  to  2  ft38  Also,  the  probes  were  placed  on  a  steel  cable  kept  under  high  tension, 


*  The  results  shown  here  differ  from  those  presented  in  Reference  [38].  First  of  all,  the 
measured  values  for  ARV-2  in  Reference  [38]  (Figures  B-4,  B-5,  and  B-6)  are  incorrect. 
Indeed,  they  appear  to  be  the  same  values  as  for  Arvonia  Shot  1,  with  the  peaks  clipped  at  2.5. 
The  computed  results  presented  here  were  also  performed  with  initial  conditions  using  the  em¬ 
pirical  constants  presented  in  Table  3-7,  while  y  =  1.3  was  used  for  all  computational  results 
presented  in  Reference  [38].  Furthermore,  the  computational  results  presented  there  for  Ar¬ 
vonia  Shots  1  and  2  used  an  initially  fiat  free  surface. 
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probe  line  was  offset  11.3  ft  from  the  charge  center  for  the  single  charge  shots.  Therefore,  we 
expect  that  the  central  plume  for  the  single  charge  shots  will  be  missed  by  the  probes,  and  only  the 
radial  plumes  will  be  detected.  The  results  for  the  single  charge  shots  BP-8,  and  BP-9,  are  shown  in 
Figures  4-23,  and  4-24,  respectively.  In  both  cases,  the  computations  indicate  larger  EWL  values  at 
earlier  times  than  the  measurements.  With  the  exception  of  the  time  interval  1.3  <  t  <  2.0  for  shot 
BP-8  in  Figure  4-23,  there  is  little  agreement  between  the  measured  and  computed  data. 

Results  for  the  multiple  charge  shot  BP-12  are  shown  in  Figure  4-25.  Since  the  computed 
values  using  both  the  P.I.  and  F.I.  quadratures  produced  similar  results,  only  the  P.I.  values  are 
shown.  The  computed  EWL  values  show  two  distinct  high  peaks  at  approximately  t  =  0.2,  and 
t  =  0.55  sec.  Peaks  at  corresponding  times  also  appear  in  the  measured  data,  but  these  are  much 
smaller  in  magnitude.  At  t  =  0.32  sec  the  top  of  the  bubble  is  actually  above  the  height  of  12.5  ft 
as  indicated  in  file  computations  shown  in  Figure  4-8.  However,  it  appears  that  the  cable  as  shown 
in  the  corresponding  photograph  has  been  significantly  deflected  upward.  This  deflection  and  rela¬ 
tively  violent  motion  of  the  cable  provides  another  possible  explanation  for  the  discrepancy  at  these 
times.  For  later  times  (t  >  0.6),  the  agreement  between  the  measured  and  computed  EWL  values  is 
very  good,  until  the  probe  data  begins  to  vanish  at  about  t  =  2.6  sec.  Indeed,  with  the  high- 
frequency  oscillations  in  the  measured  data  providing  an  approximate  bound  for  the  experimental 
error,  the  computed  data  lies  within  or  slightly  above  this  bound  for  much  of  the  interval, 
0.6  <  t  <  2.6. 

The  computed  EWL  values  shown  in  Figure  4-25,  are  in  good  agreement  with  the  computed 
microwave  data  in  Figure  4-18.  However,  there  are  important  differences  between  the  EWL  values 
determined  using  the  microwave  data  and  the  probe  data.  At  HD  =  25  ft,  Figure  4-18  shows 
EWL  >  2.5  for  0.2  <  t  <  3.2.  Since  there  will  generally  be  more  water  at  lower  elevations,  it  can  be 
expected  that  the  EWL  values  using  the  data  from  the  microwaves  at  a  height  of  12.5  ft  would  be 
substantially  above  2.5  ft  of  water,  had  they  not  saturated.  On  the  other  hand,  the  probe  data  indi¬ 
cates  that  EWL  <  1  for  a  substantial  time  in  the  interval  0.3  <  t  <  0.5.  It  therefore  appears  that  the 
magnitude  of  the  EWL  based  on  the  probe  data  is  substantially  less  than  that  indicated  by  the 
microwave  data.  Another  notable  difference  between  the  microwave  and  probe  measurements  is  the 
duration  of  the  plume.  While  the  probes  indicate  no  substantial  readings  after  3  sec,  the  microwave 
data  indicates  significant  measures  of  water  after  6  sec.  This  may  be  due  to  the  breakup  of  the 
plumes  into  small  droplets.  Since  the  probe  must  be  totally  wetted  across  the  rods,  any  droplets 
smaller  than  the  rod  spacing  will  not  be  indicated.  Furthermore,  since  the  rods  were  pointed  down¬ 
wards,  some  shadowing  could  be  expected  from  the  falling  plume.  These  limitations  are  also  dis¬ 
cussed  by  Lipton,38  who  concluded  that  the  probes  should  tend  to  underestimate  the  water  density. 

Considering  the  fact  that  the  computed  probe  data  can  be  expected  to  underestimate  the  water 
density,  the  results  for  shot  BP-14,  depicted  in  Figure  4-26,  show  significant  agreement.  The  gen¬ 
eral  profiles,  in  particular  the  time  and  duration  of  the  peak  at  time  t  ~  0.4,  and  the  overall  decay 
of  the  EWL  values,  are  in  good  agreement.  The  measured  values  are  consistently  about  50  percent 
below  the  computed  values.  The  large  computed  EWL  value  at  t  ~  0.4  is  caused  by  the  top  of  the 
water  surface  being  pushed  above  the  initial  cable  height  of  12.5  ft  This  can  be  seen  in  the 
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computed  contour  at  t  =  0.32  sec  in  Figure  4-9.  The  photograph  at  the  corresponding  time  indicates 
very  little  deflection  of  the  cable  at  this  time,  unlike  the  previous  example  where  the  agreement  was 
much  worse. 
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+  Measured  Plume  Outline 
o  Measured  Bubble  Outline 

—  Computed  p  =  0.9p0 

—  Computed  p  =  0.01po 


FIGURE  4-1.  COMPARISON  OF  COMPUTED  AND  MEASURED  BUBBLE  AND  PLUME 
OUTLINES  FOR  SHOT  CD-9  AT  TIME  t  =  0.10  SEC 
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FIGURE  4-3.  COMPUTED  AND  MEASURED  PLUME  HEIGHTS  FOR  CARDEROCK 
TEST  POND  SHOTS  CD-9,  CD-10,  AND  CD-13 
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FIGURE  4-4  COMPUTED  AND  MEASURED  PLUME  HEIGHTS  FOR  ARVONIA 
SHOTS  ARV-1  AND  ARV-2 
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(A)  INITIALLY  RECESSED  SURFACE  R{  =  1.03 d 


(B)  INITIALLY  FLAT  SURFACE 


FIGURE  4-5.  COMPUTED  PLUME  FORMATIONS  FOR  ARVONIA  SHOT  ARV-1 
(C4,  W  =  10  lbs,  d  =  6  ft)  USING  AN  INITIALLY  RECESSED 
AND  FLAT  FREE  SURFACE 
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FIGURE  4-7.  COMPUTED  AND  MEASURED  PLUME  HEIGHTS  FOR  BRIAR  POINT 
SHOTS  BP-12  AND  BP-14 
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FIGURE  4-8.  PHOTOGRAPHS  AND  COMPUTATIONS  OF  SHOT  BP-12 
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FIGURE  4-9.  PHOTOGRAPHS  AND  COMPUTATIONS  OF  SHOT  BP- 14 
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FIGURE  4-10.  EMPIRICAL  MODEL  OF  SURFACE  SHOCK  EFFECTS 
FOR  MULTIPLE  CHARGE  CONFIGURATIONS 
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FIGURE  4-11.  PHOTOGRAPHS  AND  COMPUTATIONS  OF  SHOT  BP-8 
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FIGURE  4-12.  PHOTOGRAPHS  AND  COMPUTATIONS  OF  SHOT  BP-9 
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FIGURE  4-13.  MEASURED  AND  COMPUTED  MICROWAVE  DATA  FOR  SHOT  ARV-1 
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FIGURE  4-14.  MEASURED  AND  COMPUTED  MICROWAVE  DATA  FOR  SHOT  ARV-2 
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FIGURE  4-15.  MEASURED  AND  COMPUTED  MICROWAVE  DATA  FOR  MULTIPLE 
CHARGE  SHOT  ARV-3 
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FIGURE  4-16.  COMPUTED  AND  MEASURED  MICROWAVE  DATA  FOR  SHOT  BP-8 
AT  HEIGHTS  OF  25  AND  12.5  FT 
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FIGURE  4-18.  COMPUTED  AND  MEASURED  MICROWAVE  DATA  FOR  MULTIPLE 
CHARGE  SHOT  BP-12  AT  HEIGHTS  OF  25  AND  12.5  FT 
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FIGURE  4-19.  COMPUTED  AND  MEASURED  MICROWAVE  DATA  FOR  MULTIPLE 
CHARGE  SHOT  BP-14  AT  HEIGHTS  OF  25  AND  12.5  FT 
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CHAPTER  5 

CONCLUSIONS  AND  SUMMARY 

The  primary  purpose  of  this  report  was  to  validate  the  computational  code  BUB2D  for  predic¬ 
tions  of  shallow  depth  explosion  plumes.  At  the  beginning  of  this  study  it  was  not  possible  to  reli¬ 
ably  predict  long  term  explosion  plume  behavior.  This  was  primarily  because  of  insufficiencies  in 
the  physical  modeling.  Experimental  observations  were  used  to  derive  empirical  relations,  which 
greatly  improved  the  computational  predictions.  These  modifications  to  BUB2D  are  detailed  below. 

Initial  Bubble  Conditions 

Since  the  chemistry  of  the  explosion  is  not  modeled,  empirical  laws  are  used  together 
with  the  theory  for  a  spherical  bubble  in  an  infinite  incompressible  liquid,  to  determine 
initial  conditions  for  the  BUB2D  model.  Unfortunately  all  of  the  empirical  values  for 
various  charge  types  listed  in  the  literature  were  based  on  approximations  to  the  spherical 
bubble  incompressible  theory.  These  approximations  are  valid  only  at  moderate  depths. 
New  empirical  values  are  derived  for  the  charge  types  studied  here,  without  the  use  of 
any  approximations  to  the  spherical  bubble  theory.  Furthermore,  these  values  were  based 
on  experimental  observations  at  a  variety  of  charge  depths.  Since  these  values  are  substan¬ 
tially  different  than  the  values  which  have  been  accepted  over  the  past  40  yr,  a  fairly 
detailed  justification  is  provided  in  Chapter  3. 

Bubble  Venting 

Venting  is  modeled  only  crudely  by  BUB2D,  by  instantaneously  changing  the  pressure  in 
the  vented  bubble  region  to  the  ambient  atmospheric  pressure.  It  was  found  through 
numerical  experiments  that  bubbles  would  vent  prematurely  using  the  same  cutoff  value  in 
the  code  which  has  been  successful  for  deeper  cases.  Using  a  larger  cutoff  for  liquid  cells 
near  the  water-air  interface  only,  premature  venting  was  suppresses,  while  preserving  the 
accuracy  of  predicting  the  underwater  bubble  dynamics. 

Shock  Effects 

Since  the  BUB2D  employs  an  incompressible  liquid  model,  shock  effects  cannot  be 
predicted.  The  dominant  shock  effect  for  shallow  charges  is  the  spallation  of  water  due  to 
the  shock  reflection  off  the  free  surface.  This  effect  has  been  modeled  empirically  by  ini¬ 
tializing  the  free  surface  with  an  indentation.  This  indentation  has  been  shown  to  greatly 
improve  predictions  of  the  plume  height,  particularly  at  early  times.  This  indentation  is 
characterized  by  a  single  constant  which  provided  good  plume  height  agreement  for  all 
single  charge  cases  considered. 
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The  phenomenon  of  shallow  depth  explosion  plumes  is  extremely  complicated.  The  bubble 
can  undergo  several  cycles  or  periods  before  it  vents.  During  each  cycle  hydrodynamic  instabilities 
can  be  expected  at  both  the  bubble-water  and  air-water  interfaces.  Rayleigh  Taylor  instabilities 
occur  whenever  a  denser  material  (water)  is  accelerated  into  a  less  dense  material  (bubble  gas,  or 
air).  Such  instabilities  occur  near  the  bubble  interface  when  the  bubble  is  near  its  minimum  volume 
(higher  than  ambient  pressure)  and  at  the  air-water  interface  as  the  bubble  begins  its  contraction. 
These  instabilities  appear  in  experiments  as  a  loss  of  axial  symmetry  (fingering  of  the  secondary 
plumes)  and  make  pointwise  predictions  impossible.  Nevertheless,  in  many  cases  we  were  able  to 
attain  reasonable  agreement  with  the  experiments.  A  summary  of  comparisons  of  the  computations 
to  the  measurements  are  listed  below. 

Plume  Heights 

Using  an  indented  free  surface  as  described  above,  excellent  agreement  was  obtained 
between  the  computations  and  measured  plume  heights  for  the  single  charge  cases.  How¬ 
ever,  since  we  only  used  a  two-dimensional  model  without  an  indented  surface  the  com¬ 
putations  underestimated  the  measured  plume  heights  for  the  multiple  shot  cases. 

Secondary  Plumes 

The  computed  secondary  plumes  ejected  during  the  second  bubble  pulse  are  in  good  qual¬ 
itative  agreement  with  the  observations.  For  many  computational  models  this  is  a  difficult 
feature  to  compute  since  it  occurs  well  after  the  first  bubble  collapse.  In  addition  to  the 
relatively  long  time  it  takes  to  develop,  the  phenomenon  is  further  complicated  by  the  fact 
that  the  bubble  may  form  one  or  more  toroidal  regions  after  this  collapse.  Neither  of 
these  factors  is  an  issue  for  the  model  in  BUB2D  and  no  special  modifications  to  the  code 
were  made  concerning  the  prediction  of  these  plumes. 

Microwave  Data 

Considering  the  complexity  of  the  problem  and  the  lack  of  an  independent  validation  for 
the  measurements,  comparisons  between  the  computations  and  experiments  were  reason¬ 
ably  good.  Despite  the  fact  that  the  multiple  shot  plume  heights  were  underpredicted, 
there  was  remarkable  agreement  between  the  computed  and  measured  rise  in  EWL  values 
at  early  times.  This  suggests  that  the  visible  peaks  of  the  plumes  contain  a  relatively 
insignificant  mass  of  water. 

Density  Probe  Data 

The  only  significant  match  between  the  measured  and  computed  probe  data  occurred  with 
the  multiple  shot  cases  BP-12  and  BP-14.  The  probes  must  be  sufficiently  close  together 
and  the  data  must  be  integrated  to  provide  any  meaningful  comparisons. 

One  issue  not  addressed  in  detail  in  this  report  is  the  three-dimensional  effects  that  occur  when 
using  multiple  discrete  charges.  Several  experiments  were  conducted  at  Dynamic  Testing  Incor¬ 
porated  in  Rustberg,  Virginia  in  June,  1995  to  investigate  these  effects.  This  data  will  be  used 
together  with  BUB3D  and  the  empirical  model  presented  in  Figure  4-10  in  a  forthcoming  publica¬ 
tion. 
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The  algorithm  used  for  both  codes  BUB2D  and  BUB3D  are  continually  being  improved.  One 
such  improvement  is  to  ensure  the  strict  conservation  of  mass.  A  version  of  the  codes  that  do  not 
add  mass  when  “numerical  cavitation”  occurs  (see  Chapter  2)  is  currently  being  tested.  Since  the 
mass  is  the  most  important  variable  for  these  shallow  plume  predictions,  this  modification  can  be 
expected  to  improve  our  predictive  capability. 

Based  on  the  validations  presented  in  this  study,  it  is  now  possible  to  conduct  optimality  stu¬ 
dies  with  some  degree  of  confidence.  For  example,  it  may  be  possible  to  optimize  the  duration  that 
the  EWL  is  above  some  prescribed  threshold  value.  Such  studies  will  be  the  topic  of  future  work. 
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